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Summary 


Measurements of the heat transfer, recovery factor, and 
pressure distributions around a circular cylinder normal to a 
supersonic rarefied-air stream (total temperature = 300°K.) 
are described for the Mach number range of 1.3 to 5.7, the 
Reynolds number range of 37 to 4,100 and at cylinder wall 
average temperature levels of 90°K. and 210°K. Study of the 
results yielded: (1) a correlation equation for the stagnation- 
point Nusselt number as a function of the Reynolds number just 
after the normal part of the detached bow shock wave; and (2) 
Fourier series expressions for the heat-transfer coefficient and 
pressure coefficient distributions in terms of the stagnation point 
values. 

In comparing these measurements with predictions based on 
recent analytical studies, exceptionally good agreement for the 
heat-transfer coefficient distribution was obtained with Lees’ 
theory.!! In the Mach number range of 3.55 to 5.73 the pres- 
sure decreased less rapidly with distance from the stagnation point 
than predicted by the modified Newtonian theory. 


Symbols 


a,a’,c = Fourier coefficients 


A = area (ft.?) 

Sy = specific heat at constant pressure (B.t.u./Ib. °F.) 

C,(@) = coefficient of pressure at angle 6; 
= (Pw — po)/(1/2)pata? 

C,(0) = coefficient of pressure at the forward stagnation 
point 

d = cylinder diameter (ft.) 

h(@) = heat transfer coefficient at angle @(B.t.u./hr. ft.2 °F.) 


Received July 29, 1959. Also presented at the Hypersonic 
Phenomena Session, IAS 28th Annual Meeting, N.Y., Jan. 25-27, 
1960. 
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Heat Transfer, Recovery Factor, and Pressure 
Distributions Around a Circular Cylinder 
Normal to a Supersonic Rarefied-Air Stream‘ 


O. K. TEWFIK* anp W. H. GIEDT** 
Unwersity of California at Berkeley 


h(O) = heat transfer coefficient at the forward stagnation 
point (B.t.u./hr. ft.2 °F.) 

ks = thermal conductivity of material of cylinder 
(B.t.u./hr. ft. °F.) 

de = heat transfer rate by convection (B.t.u./hr.) 

qrad = heat transfer rate by radiation (B.t.u./hr. ) 

r = radius (ft.), distance from origin in polar coordi- 
nates 

R = recovery factor = (Taw — To)/(T1.. — To) 

t = temperature (°F.) 

bey = wall temperature (°F.) 

Taw = adiabatic wall temperature (°F. ) 

T 2 = free stream static temperature (°K.) 

Fen = free stream total temperature (°K.) 

u = stream velocity (ft./sec. ) 

Ue = tangential velocity component at outer edge of 


boundary layer (ft./sec. ) 


0 = angle measured from stagnation point streamline 
(degrees ) 

p = density (Ib./ft.*) 

= viscosity (Ib./sec. ft.) 

Ue = tangential velocity component at outer edge of 
boundary layer, in transformed coordinates 

X = tangential transformed Cartesian coordinate 

Subscripts 

oo = free stream 

1 = just behind normal shock wave 

e = at outer edge of boundary layer 

w = outer surface of test model 


(I) Introduction 


Sng FLOW FIELD around a circular cylinder normal 
to a supersonic air stream is characterized by the 
presence of a detached bow shock wave ahead of the 
cylinder and a boundary layer around most of its 
circumference. Due to the viscous nature of the flow 
in the shock wave and the boundary layer, part—and 
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Fic. 1. Thick-walled heat transfer cylinder. 


in some places all—of the kinetic energy of the air 
stream is converted into enthalpy. The result is a 
high temperature region around the cylinder, par- 
ticularly in the neighborhood of the forward stagnation 
point, from which heat may be convected and radiated 
to the surface. When the free stream temperature is 
less than about 400°K. and its Mach number less than 
6, no dissociation of the gas molecules will occur in or 
behind the bow shock wave, and radiation from the 
gases will be negligible. If, in addition, the free 
stream is rarefied, its Reynolds number will be low 
enough in most cases for the boundary layer to be 
laminar over the cylinder surface up to the separation 
point. Under these conditions, a theoretical speci- 
fication of the heat transfer distribution from the 
forward stagnation point up to the separation point 
involves the study of compressible laminar boundary 
layer theory with streamwise pressure variation. 

To obtain an exact solution to this problem is a 
formidable undertaking, and has proved impossible 
so far. Consequently, various different assumptions 
and simplifications have been introduced to obtain 
solutions, the most common being one or more of the 
following: (1) Pr = 1, (2) Pr = constant, (3) the heat 
transfer surface is isothermal, (4) the velocity at the 
outer edge of the boundary layer varies as a simple 
power function of the coordinate parallel to the flow, 
(5) the concept of local similarity, (6) the concept of 
a “highly cooled surface,’’ and (7) the viscosity varies 
as the first power of the temperature. Solutions to the 
problem thus simplified have been reported for the 
laminar compressible boundary layer equations with 
heat transfer and pressure variation, and may be 
applied to cylinders in cross flow as well as to axially 
symmetric bodies. 

Comparatively few experimental investigations to 
test these solutions have been reported, and most of 
these have mployed hemisphere-cylinder models 
with values of© Re. (based on model diameter) on the 
order of 104-1 ° and values of M., in the range 1.5—4 
(see, e.g., references 1 and 2). Therefore, to augment 
the meager amount of experimental data available, and 
to extend the measurements into the lower ranges of 
Re.,, an experimental investigation of the local heat 
transfer around a cylinder normal to a rarefied-air 
stream was recently conducted in the University of 
California low pressure wind tunnel. A circular cylin- 
der normal to the stream flow was selected for the 


model in view of the following: (1) local heat transfer 
to such bodies is of both fundamental and technical 
importance; (2) conditions on the forward half are 
characteristic of those on a blunt-nosed body, a subject 
of current interest; and (3) local heat transfer meas- 
urements are inherently difficult and a cylindrical model 
offers some advantages with regard to necessary in- 
strumentation. The experimental method is described 
in the following section. The test results are then pre- 
sented and compared with several pertinent theories. 


(II) Method and Scope 


With interest centered on aerodynamic heating it was 
desired to have heat transfer from the tunnel air stream 
to the model. To provide this condition it is necessary 
either to heat the air stream, or to cool the model. The 
latter method was adopted, because the available 
tunnel facility was not designed to operate under high 
temperature conditions. The cylindrical model was 
cooled to average temperatures around 90° and 210°K. 
by pumping liquid nitrogen and butyl cellulsolve 
(cooled by dry ice), respectively, through it. With a 
tunnel free-stream total temperature around 300°K.., 
the resulting wall-to-stream temperature ratios were 
0.3 and 0.7, which correspond to flight at Mach num- 
bers of 4 and 1.5. 

High-speed convective heat transfer is conventionally 
specified in terms of Newton’s law of cooling, except 
that the adiabatic wall temperature is used in place 
of the fluid stream static temperature; that is, 


(q-/A) = K(taw tw) (1) 


Under steady-state conditions heat is transferred to the 
outer cylinder surface by convection and radiation, the 
sum of which must be equal to the heat conducted 
inward through the cylinder wall. Thus, 
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Fic. 2. Pressure coefficient distributions around adiabatic and 
cooled cylinder for M. = 5.73. 
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(qe/A) + (Graa/A) = — [k.(Ot/Or)],, (2) 
Substituting from Eq. (1) and solving for h gives 


tow — ty 

Equation (3) indicates the basis of the experimental 
approach used in this investigation. To compute the 
heat conduction through the cylindrical model at 
various circumferential locations, the inner and outer 
surface temperature distributions were determined 
from resistance measurements of thin metallic films* on 
the inner and outer surfaces of the electrically noncon- 
ducting (glass) model (Fig. 1). These films were 
located at the same angular position and were only 
8° of are wide, thus giving essentially local tempera- 
tures. The circumferential temperature distributions 
were determined by rotating the cylinder. Circum- 
ferential heat conduction was accounted for by evalu- 
ating (Of/Or),, from a numerical series solution of the 
Laplace equation for the temperature distribution in 
the cylinder wall. This solution was determined from 
the temperature measurements on the inner and outer 
surfaces taken at small angular intervals by rotating 
the cylinder in the airstream. Values of k, for the glass 
used were available. (g,aa/A) was evaluated from 
the temperature differences observed with no airflow 
across the cylinder. 

A special thin-walled cylinder was constructed for 
measuring the adiabatic wall temperature, f,,. Its 
wall thickness was made about 0.005 in. to reduce 
circumferential conduction to a minimum. Tempera- 
tures were determined from resistance measurements 
of a thin metallic film on its outer surface. 

Pressure measurements were made with an 0.5-in. 
O.D. copper tube which had a pressure tap 1 in. long 
and 0.010 in. wide milled along a generator. The 
slot was connected to the pressure sensing device 
through an 0.25-in. O.D. tube inside the 0.5-in. O.D. 
tube. Thus liquid nitrogen could be pumped through 
the annular space to cool down the model during pres- 
sure distribution measurements, if so desired. 

A sensitive Wheatstone bridge was specially built 
to measure the resistance of the metal film strips. It 
could easily detect changes of resistance of 1 part in 
100,000. 

For more details regarding the development of the 
heat transfer model and the calibration of its film strips, 
the construction of the adiabatic wall temperature 
model, the pressure distribution model, and the model 
rotating mechanism, reference 5 may be consulted. 
Test procedure, data reduction, and results in tabular 
form are also reported there in detail. 


(III) General Characteristics of Results 


(A) Pressure Distribution 


Representative pressure coefficient distributions for 
both cooled and uncooled models are plotted in Figs. 
2 and 3. In general, the data defined very smooth 
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Fic. 3. Pressure coefficient distributions around adiabatic and 
cooled cylinder for M. = 1.32. 


curves. There is very little difference (about 5 per cent 
or less, attributable to experimental scatter) between 
the pressure distributions for the cooled and adiabatic 
conditions, except in the wake region at low values of 
Re. (see Fig. 3). 


(B) The Boundary Layer Separation Point 


The location of boundary layer separation from a 
cylinder surface can usually be determined by the pres- 
ence of an inflection point in the pressure distribution 
curves. In the present study, however, no inflection 
point was evident in the C,(@) curves, despite all the 
experimental care taken to look for it. The boundary 
layer was relatively thick, due to comparatively low 
values of Re.., and so could not be expected to remain 
attached very long in an increasing pressure region. It 
was therefore decided to associate separation with the 
angle 8, which denotes the location of the minimum 
value of C,(@). In general, 8 varied from 120° to about 
160°, which is appreciably higher than in subsonic flow. 
Qualitatively, 8 increases as Re., decreases, as M, 
increases. For the same values of the free stream 
parameters, 8 is always larger for an adiabatic as com- 
pared to a cooled surface. 


(C) A General (Empirical ) Formula for C,(@)/C,(0) 


It was noted that the curves of C,(@) for all the dif- 
ferent free stream Mach numbers had the same general 
shape, and were almost identical at the higher values 
of M.. This suggested the possibility of describing 
the pressure distribution in terms of one general equa- 
tion valid for the entire cylinder surface. From sym- 
metry considerations, this general expression should be 
an even function of 6. Since a particular location on 
the cylinder surface may be designated by any of the 
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angles 0,0 + 27,0-+ where 
0 < 6 < 7, the general expression has to be a periodic 
function of 6 with period 27. The only function of 6 
which satisfies these two conditions is a Fourier cosine 
series. Hence it should be possible to describe C,(@) 
for each test condition by a series of the form do’ + a,’ 
cos 6 + cos 2(0) + a;’ cos 3(0) + ...forO 
Determining the coefficients a,’ showed that a,’ was 
less than 0.002 a,’ in all cases; this term was therefore 
neglected, along with subsequent terms. The resulting 
expression for C,(@) was thus 


C,(0) = ao’ + ay’ cos 6 + a,’ cos 20 + a;’ cos 36 (4) 


To compare the results for the various free stream 
conditions, the ratio C,(@)/C,(0) was determined. 
Since from Eq. (4) C,(0) = (ado’ + ai’ + a2’ + a;’), 


C,(0)/C,(0) = ao + a: cos (6) + 
az cos 26 + a3 cos 36 (5) 


where dp = do’/(ao’ + ai’ + a2’ + a;’), etc. For 


T T T T T T T T T T T U T T 
—— EMPIRICAL FORMULA FOR HYPERSONIC FLOW 
Cpl@) 
= 0.30 + 0.46 cos 8+0.20c0s28 +0.04 cos38 
Cp(0) 
| © EXPERIMENTAL,ADIABATIC MODEL, Mq@=5.73 
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Fic. 4. Comparison of pressure distribution predicted from 
empirical equation with data points for 7. = 5.73. 
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Fic. 5. Typical recovery factor distributions. 
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Fic. 6. Heat transfer correlation at stagnation point. 


> 4, the coefficients a, d2, become more or 
less independent of /.,.. Values of a, varied by less 
than 3 per cent from their mean values. Substituting 
these mean values in Eq. (5) gave 


C,(0)/C,(0) = 0.30 + 0.46 cos 0 + 
0.20 cos 29 + 0.04 cos 36 (6) 


which is accurate to within 3 per cent for all values of 
M.. => 3.5, and independent of heat transfer to the 
model. Comparison with the data from a typical 
experimental run (/7., = 5.73, Re. = 4,100) is shown 
in Fig. 4. For values of 7. < 3.5, no obvious trend 
was apparent in do, @1, dz, and 


(D) Recovery Factor 


Typical recovery factor distributions are plotted in 
Fig. 5. At the stagnation point the mean value of all 
the measurements is 1, and the standard deviation is 
+0.002; at @ = 45°, R has a mean value of 0.973, with 
a standard deviation of +0.006 (the primary reason 
for the experimental scatter was the slight variation in 
total free stream temperature during the course of each 
run). Thus, for 0 < 6 < 45°, the recovery factor may 
be considered independent of free stream Mach or 
Reynolds numbers. Beyond 6 = 45°, R in general 
decreased to a minimum of about 0.91 around 135° 
and then increased slightly to about 0.93 at the rear 
stagnation point. Investigation of the effect of radia- 
tion from the tunnel walls showed, however, that the 
value of R in the rear stagnation region could be high 
by as much as 25 per cent at the lower values of Re...° 


(E) Heat Transfer at the Forward Stagnation Point 


In considering the possible correlation of the heat 
transfer coefficient at the stagnation point in terms of a 
Nusselt number as a function of Reynolds and Mach 
numbers, some effect of the two different temperature 
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Fic. 7. Heat transfer correlation at stagnation point for hyper- 
sonic flow. 


levels (average surface temperatures of 90° and 210°K.) 
employed was anticipated. Various procedures were 
therefore investigated. The first involved computing 
the Nusselt number with the thermal conductivity of 
the air stream evaluated at the wall stagnation point 
temperature, and the Reynolds number with the 
viscosity, density, and velocity for conditions just 
behind the normal shock. When plotted, these values 
fell along two distinct curves, one for each coolant. 
The next step was to compute Reynolds number with 
the viscosity and density based on wall temperature. 
Again the plotted results indicated two curves, with 
considerable scatter around each. The third and most 
successful combination investigated was to base both 
Nusselt and Reynolds numbers .on conditions just be- 
hind the normal shock, i.e., to define the numbers as 
follows: Nu,(0) = hd/ki, Rey = d mpi/ms. Nu(0) 
is plotted versus Ke; in Fig. 6 for all free stream con- 
ditions and the two levels of wall temperature. The 
expression 


Nu,(0) = 0.38 Re,°-* (7) 


represents the data points within +6 per cent. The 
fact that the results for the two temperature levels 
could be correlated without introducing some effect 
of the different temperature levels means that the heat 
transfer throughout the temperature range involved 
is essentially independent of fluid property variation. 
From Eq. (7) it is apparent that the stagnation point 
Nusselt number in supersonic flow is a function of the 
Reynolds number, with all properties evaluated at 
conditions just behind the normal part of the detached 
shock wave ahead of the cylinder. The influence of 
M.. is not directly manifested, although it determines 
the strength of the normal shock. This suggests an 
alternate relation may exist for Nw,(0) in terms of 


Re, and M.,. To determine such an expression it was 
assumed that 7)/7., ~ M..*, which is approximately 
correct for hypersonic values of M/,,. Also assuming 
that ~ 


Re:/Re. ~ 1/M.'* (8) 


Substituting Re, from Eq. (8) in Eq. (7) and rearranging 
yielded 


Nu,(0)-M., ~ (9) 


To verify the above equation, and to evaluate the con- 
stant of proportionality, values of Nw,(0)\/.. were 
plotted versus Ke,, for the hypersonic range (Fig. 7). 
The results are described by the equation 


Nu(0)-M., = 0.67 (10) 


within +7 percent. This agreement supports the orig- 
inal correlation, Eq. (7), and the validity of Eq. (9). 


(F) Distribution of Local Heat Transfer Coefficient 


Although it is possible to compute values for local 
Nusselt numbers at the angle @ at which values of the 
local heat transfer coefficient /4(@) were obtained, the 
length and the value of the thermal conductivity on 
which they are based are somewhat arbitrary, and as a 
result, of dubious significance. For this reason it was 
decided to consider the heat transfer coefficient dis- 
tribution directly. 

Typical results are shown in Fig. 8. As can be 
noted maximum values occur at the forward stagnation 
point with a gradual decrease to almost negligible 
values on the rear half of the cylinder. In view of the 
success in describing the pressure distributions by a 
Fourier series, a similar procedure was applied to the heat 
transfer data. That is, 4(@)/h(0) was approximated 
by a series of the form: 
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|— n(6)/n(0) =0.37+.0.48 +0.15c0s 28 
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Fic. 9. Comparison of local heat transfer distribution at 
M. = 5.50 with empirical equation. 


h(6)/h(O) = c + cos + 


cos 20 + c;cos 306+... (OS (11) 


Values of c; proved to be about 2 per cent or less of &, 
in all except two cases (where c; was 20 and 25 per cent 
of c,); consequently only the first three terms of the 
series were retained. The respective coefficients were 
nearly equal to each other and, taking the mean, the 
expression 


h(0)/h(0) = 0.37 + 0.48 cos 6 + 


0.15 cos 20 (0O<@< 7m) (12) 


was found to represent the experimental data within 
+6 per cent. A graphical comparison with actual 
data points for /., = 5.5 is shown in Fig. 9. 

Since this expression is based upon 20 different tests 
with widely differing free stream conditions and two 
significantly different temperature levels, it may be 
concluded that 4(@)/h(0) is essentially independent of 
free stream conditions or wall temperature level 
throughout the range of this investigation. 


(IV) Comparison of Experiment With Theory 


(A) Pressure Distribution 


The pressure distribution predicted from the modi- 
fied Newtonian theory® is compared in Fig. 10 with the 
empirical formula which was found to best fit the ex- 
perimental pressure distributions. As may be seen, 
the experimental curve does not decrease with dis- 
tance from the stagnation point as rapidly as the theo- 
retical. (The theory does not predict any values for 
68> 90, so the comparison has to be limited to0 < 6 < 
90, although measurements were continued around the 
entire cylinder surface.) 

In seeking an explanation for this deviation it is 
noted that the values of Re, in the present tests were 
substantially smaller than in previous tests! in which 
rather good agreement with the modified Newtonian 
theory was found. It is likely that the thicker bound- 
ary layers associated with these lower values of Re. 
caused the effective shape of the test model to differ 
significantly from a circular cylinder. This hypothesis 
is further supported by the fact that the boundary 
layer growth on axially symmetric blunt-nosed bodies 


(where good agreement with the Newtonian theory is 
obtained, even at low values’ of Re.) is much slower 
than on two-dimensional blunt-nosed bodies when 
submerged in the same free stream, and therefore 
has less effect on the pressure distributions. 


(B) Heat Transfer at the Forward Stagnation Point 


The recent analytical investigations of heat transfer 
from a laminar compressible boundary layer by Stine 
and Wanlass,* Cohen and Reshotko,® and Lees!! 
were considered particularly applicable to the present 
experimental study. In making a comparison it was 
found advantageous to consider first the stagnation 
point and then the variation of the heat transfer 
coefficient relative to its stagnation point value. 

Cohen and Reshotko present their results” for the 
stagnation point in terms of the heat transfer parameter 


_ h(0)x/Rew 
V Rey V x ptte/ 


as a function of 7,,/7;,,. with Pr, asa parameter. For 
each of the twenty heat transfer tests conducted, 
T/T... and Pr, were computed and Nu,,(0)/ V Ren 
was evaluated from their curves. It was then inter- 
esting to note that although 7,,/7;,. varied from about 
0.24 up to 0.75 in the experiments, Nu,,(0)/ V Rey 
as predicted from the theory was practically constant 
and equal to 0.48. The reason is that when 7, is in- 
creased, Pr, is decreased so that Nu,,(0)/ V Rew 
remains almost unchanged. 

It is preferable to define Nu, and Re, for the experi- 
mental results in terms of d rather than x, since both 
parameters become zero at the stagnation point if x 
is used, making the ratio Nu,,(0)/ V Re, indeterminate. 
If it is assumed that u,(x) = Cx near the stagnation 
point, where C is a constant, 
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Fic. 10. Comparison of empirical pressure distribution equation 
for hypersonic flow with the modified Newtonian theory. 
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Nu,,(0) h(0)d/Ry 

(13) 

where i,’(0) = dC/m is the dimensionless velocity 


gradient at the stagnation point. Equating this quan- 
tity to 0.48 leads to the following expression for com- 
parison with the experimental data: 


Nutya(0)/W = 0.48V Rew.a (14) 


The experimental values of Nuy.¢ (0)/ V ii,’(0) are 
plotted versus Re,,, in Fig. 11, together with the theo- 
retical relation (14). The agreement is fair, the stand- 
ard deviation of the experimental results from the theo- 
retical line being about 14 percent. 

The experimental results may also be compared 
with Lees’ approximate expression for the heat transfer 
at the stagnation point,'! 


Nu,(0) = h(0)d/k, = 0.50 Pr’? V (15) 


which is based on the concept of a highly cooled surface. 
Assuming that u ~ 7, p ~ k ~ T, (these condi- 
tions are very closely satisfied at the low temperatures 
of the experiments), this expression becomes inde- 
pendent of the temperature. Taking Pr = 0.7, it 
finally reduces to: 


Nu(0)/Wi,'(0) = 0.44V Re, (16) 


This is plotted in Fig. 12 and compared with the experi- 
mental results obtained when the model was cooled 
with liquid nitrogen (model surface temperature around 
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Fic. 11. Stagnation point heat transfer; comparison of Cohen 
and Reshotko’s theory with experimental data points. 
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Fic. 12. Stagnation point heat transfer; comparison of Lees’ 
theory with experimental data points. 


90°K., and p»/p, ~ 3) and with butyl cellusolve cooled 
with dry ice (model surface temperature around 
210°K., and p,y/p, ™ 1.5). The standard deviation 
of the data from the theoretical line is about 9 per cent, 
which is a substantial improvement over that found with 
the curve due to Cohen and Reshotko. 

It should be noted that Lees’ theory assumes that 
Pu/p. be much larger than 1, say 10-20. This condi- 
tion was introduced for the case of hypersonic Mach 
numbers, where it is very well satisfied. In the present 
experiments, however, p,/p- was 3 in one series of tests, 
and only 1.5 in the other series, and thus much lower 
than postulated. It is therefore remarkable that the 
good agreement shown in Fig. 12 was obtained. 

It was observed that Lees’ equation is very similar 
to Squires’ formula,'? derived on the assumption of 
constant fluid properties: 


Nu(0) = 0.570 (Pr)* VW Re, (17) 


where i,’(0) = 4 for ideal fluid flow around circular 
cylinders. Substituting the value 0.70 for Pr, Eq. (17) 
becomes: 


Nux(0)V = 0.495V Re: (18) 


which is about 10 per cent higher than Eq. (16). Thus, 
if only the numerical factor in the above equation is 
reduced by 10 per cent, the resulting formula will fit 
the experimental results within 9 per cent. This is, 
of course, consistent with previous findings regarding 
the negligible effect of fluid property variation with 
temperature. 

In comparing the theoretical Eqs. (14) and (16) with 
the empirical Eq. (7), it is interesting to note that the 
Nusselt number at the stagnation point varies as Re,°° 
in the former and as Re*- in the latter. This apparent 
anomaly may be resolved by noticing that in Eqs. (14) 
and (16) Nu(0) also depends on #,’(0), which was not 
constant in the ranges of M., and Re. investigated, due 
primarily to the following: 

(1) Inviscid flow effects: As M,, increases from 1.32 
to 5.7, the pressure gradients at the outer edge of the 
boundary layer become larger with a corresponding 
increase in i,’(0). Since Re, increases with M,, in 
the wind tunnel used, #,’(0) increases with Re... 
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Fic. 13. Heat transfer coefficient distribution; comparison 
theory with empirical equation. 


(2) Viscous flow effects: As Re. increases from 37 
to 4,100, the boundary layer becomes thinner, and its 
interaction with the free stream less pronounced. The 
“effective body’’ becomes less blunt, and 7#,’(0) in- 
creases. 

As a result of the combined effects of these two fac- 
tors, 7,’(0) increased from 2.5 to 5.4 when Re, changed 
from 37 to 380 as computed from the measurements. 
Thus i,’(0) ~ Re,'/*, and when this is substituted in 
Eqs. (14) and (16), the result is that Ni(0) ~ Re,-*, 
which then agrees with the correlation Eq. (7). 


(C) Heat Transfer Distribution 


This may be specified by conventional laminar com- 
pressible boundary layer theory with streamwise pres- 
sure variation. Considering first the work of Stine 
and Wanlass,* it is noted that they assumed linear 
variation of viscosity with temperature, constant c, 
and Pr, and that they used the Stewartson transforma- 
tion to transform the momentum and energy equa- 
tions to the incompressible form. In the transformed 
coordinates the tangential velocity at the outer edge 
of the boundary layer U’, was assumed proportional to 
the distance from stagnation point X according to a 
power law, i.e., U.~X"™, where m isa constant. The 
transformed momentum and energy equations were 
then uncoupled by the assumption of a small tempera- 
ture potential compared with the free stream total 
temperature. Solutions of the two equations for 
various values of the index m are known (see, for ex- 
ample, reference 13). 

Since in the case of a cylinder in cross flow U, does 
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not vary as X”™, the concept of local similarity was 
introduced to specify the local temperature gradient. 
The idea involved may be stated as follows: in the 
transformed coordinates, the dimensionless tempera- 
ture gradient at the cylinder surface at a particular 
point X where the velocity and the velocity gradient 
at the outer edge of the boundary layer are U, and 
dU,/dX, is identical with that at the same value of X 
on the surface of a wedge, the flow over which is such 
that the distribution of the velocity U, at the outer 
edge of the boundary layer is given by UL’, ~ X™ where 
m is a constant equal to (X/U,)(0U,/dX). Thus only 
the local values U, and dU,/dX at a particular point 
X are used to determine the heat transfer there, irre- 
spective of the previous variation of LU’, with X up to the 
point under consideration. 

From the theory, 4(@)/h(0) was computed for the 
special case of /., = 1.90, using the measured pressure 
distribution. The results are compared in Fig. 13 
(curve 1) with the empirical Eq. (12) which describes 
the results for all free stream conditions and two differ- 
ent wall temperature levels. The theoretical curve 
predicts a much slower decrease from the stagnation 
point than was found experimentally, the deviation 
being approximated by 0.4{1 — 
believed the reason for this is that the boundary layer 
development on the cylinder is appreciably different 
from that determined by the wedge flows associated 
with the various locations on the cylinder surface. 

The analysis of Cohen and Reshotko® is very similar 
to that of Stine and Wanlass except that they solved 
the transformed energy and momentum equations 
simultaneously instead of uncoupling them by the as- 
sumption of a small temperature potential. They as- 
sumed Pr = 1, and accounted for viscosity variation 
with temperature with the Chapman-Rubesin coeffi- 
cient A, ie., w/uo = AT/T>. Their predictions thus 
show slightly better agreement (curve 2 in Fig. 13) 
with experiment than Stine and Wanlass found, the 
deviation from the empirical equation being about 0.3 
{1 — [h(0)/h(0)].2p}. It may be inferred in passing 
that the assumption of small temperature potential 
does not appreciably affect the heat transfer distribu- 
tion, and hence the momentum and energy equations 
may be uncoupled and their solution greatly simplified 
without introducing substantial error. 

Lees'' used Levy’s transformation to reduce the 
laminar compressible boundary layer equations to the 
incompressible form. Then the concept of the highly 
cooled surface was introduced, (p,/p» < 1), such as 
applies in flight at hypersonic speeds, when the surface 
is at a very low temperature compared with the free 
stream total temperature, to prevent it from disin- 
tegrating. With this assumption Lees reasoned that 
the pressure gradient term in the transformed momen- 
tum equation has negligibly small affect on the local 
dimensionless enthalpy gradient, and may therefore 
be omitted. His methods were followed in specifying 
the distribution of (@)/h(0) for W., = 1.90 and 4.19. 
The results, which are plotted as curves 3 and 4 in 
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HEAT TRANSFER, 


Fig. 15, showed rather good agreement with the experi- 
mental distribution, the deviations from the experi- 
mental being about +0.1}1 — [h(6)/h(0)| 

That such good agreement with the experimental re- 
sults was obtained following the method of Lees sug- 
gests that the influence of the pressure gradient is minor 
even when the condition p,/p» << 1 is not satisfied. In 
view of this it may be possible to achieve some im- 
provement in the theory by simply selecting a suitable 
constant value for the pressure gradient term, as in 
stagnation point flow, rather than neglecting it com- 
pletely. 


(V) Summary and Conclusions 


(1) The pressure coefficient distribution, when related 
to its stagnation point value, becomes independent of 
the free stream Mach and Reynolds numbers or heat 
transfer to the model when the free stream Mach num- 
ber exceeds about 3.5, and is given within + 3 per cent 
by the empirical equation: 


C,(@)/C,(0) = 0.80 + 0.46 cos + 
0.20 cos 20 + 0.04 cos 36 


The above expression, when compared with the modi- 
fied Newtonian theory, is larger by about 0.1}1 — 
[C,(0)/C,(0)] theor } 

(2) The recovery factor at the forward stagnation 
point is 1.0 and decreases slowly along the cylinder 
circumference to a minimum value of about 0.91 
around @ = 135° and then rises slightly as @ approaches 
180°. It is essentially independent of free stream 
Mach and Reynolds numbers. 

(3) The Nusselt number at the forward stagnation 
point is best correlated with the Reynolds number based 
on conditions just behind the normal shock; the equa- 
tion 


Nu,(0) = 0.38 


represents the data within +6 per cent. In comparing 
the data with analytical predictions, the best agree- 
ment was with Lees’ theory, being +9 per cent. 

(4) The heat transfer coefficient distribution in 
terms of the stagnation point value is independent of 
free stream Mach and Reynolds numbers, and wall 
temperature level in the ranges investigated, and is 
given within +6 per cent by the empirical equation: 


RECOVERY FACTOR 729 


h(0)/h(0O) = 0.37 + 0.48 cos 6 + 0. 5 cos 26 


When compared with the predictions from three par- 
ticularly applicable recent theories, the best agreement 
was obtained using Lees’ approach, in spite of the fact 
that the cylinder surface was not “highly cooled’’ as 
assumed in the analysis. This assumption is appar- 
ently not necessary; it may be relaxed altogether with- 
out affecting the final results or introducing any addi- 
tional complications. 
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Second-Order Theory for Unsteady 


Supersonic Flow Past Slender, 


Pointed Bodies of Revolution‘ 


JAMES D. REVELL* 


North American Aviation, Inc. 


(1) Summary 


An analysis is made of the second-order effects of thickness on 
the unsteady aerodynamic forces on a slender pointed body of 
revolution in supersonic flow. The theory is restricted to har- 
monic oscillations for small angles of attack. The solution is ob- 
tained by approximating the nonlinear terms in the second-order 
potential equation by their first-order values and solving the 
resulting inhomogeneous partial differential equation, subject to 
more refined boundary conditions. The pressure equation is 
likewise refined and integrated to give the second-order correc- 
tions to lift and pitching moment coefficients. The analysis can 
be considered as an extension of the second-order, slender body 
theory of Lighthill to the case of unsteady flow. 

The results indicate appreciable reductions in unsteady lift 
and damping moment coefficients when applied to slender cones. 
The present theory is estimated to be reliable provided that M6 
is less than 0.7. 


(2) Symbols 


U = flight velocity, ft./sec. 
= dimensionless local acoustic velocity = 1/M/ 
l = body length, ft. 
Q = dimensionless total velocity potential 
T = dimensionless time measured in a Newtonian 


reference frame 
= specific heat ratio [also Eulers constant, 0.5772... 


Y 
in Section (5)] 

a = instantaneous angle of attack due to harmonic 
pitching 

= dimensionless body axis coordinates 

1,,1,, 16 = unit vectors in cylindrical body coordinate sys- 
tem 

6 = 2e = dimensionless base diameter 

t = dimensionless time measured in body coordinate 
system 

v,? = two-dimensional Laplacian operator (dimension- 
less) 

M = free-stream Mach number 


o = velocity perturbation potential (dimensionless ) 
¢ = first-order slender body potential 

Cy = second-order slender body potential 

= second-order slender body correction potential 
y = particular solution to V 

x = complementary solution to v 

V, = dimensionless body velocity 
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G = pressure coefficient, C, = (p — p.)/(1/2)p. U? 

Cy, Cn = coefficients of Fourier series for C, 

B = 

Po = static pressure of undisturbed fluid, lbs. per sq.ft. 

Bw = density of undisturbed fluid, slugs/cu-ft. 

R(x) = dimensionless body profile radius 

S(1) = dimensionless base area = ze? 

F = dimensionless pressure function, F = 20, + 
(VQ)? — 1 

w = angular frequency, rad./sec. 

k = reduced frequency, k = wi/U 

XG = axis of rotation (dimensionless) 

Xo = reference axis for moments (dimensionless ) 

Ko, Ki = modified Bessel functions of the second kind of 
order zero and one 

s = Laplace transform complex variable 

q; = modified Bessel function of the first kind of order 
one 

L( ) = Laplace transform of function with respect to x 

K = hypersonic similarity parameter 


Subscripts and Superscripts 


0 = axial flow component 

1 = cross flow component 

ro) refers to the properties of the undisturbed air 

H refers to harmonic variation of a; also 
homogeneous 

VA refers to constant vertical acceleration 

( )® = remainder 

c~) = Laplace transform of quantity 

t..) = vector quantity 

( subscripts denote partial derivatives (except 


when used with unit vectors) 
C¥ primed variables are dimensional values in 
connection with velocities and coordinates 


(3) Introduction 


hase to the first-order potential equation for a 
harmonically oscillating, slender, pointed body of 
revolution have been obtained by Miles! and Dorrance.’ 
Sacks’ has extended Miles’ solution and has evaluated 
stability derivatives for bodies of general cross-sectional 
shapes. Miles’ and Sacks’ methods are an extension 
of the well-known slender body theory of Ward‘ to the 
unsteady case. The solution of Dorrance*® is based 
on the source integral equation approach and utilizes 
slender body concepts to satisfy the boundary condi- 
tions on the body. 

The purpose of the present paper is to investigate 
the second-order effects of thickness on the unsteady 
supersonic flow past slender pointed bodies of revolu- 
tion by means of an iteration procedure wherein the 
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nonlinear terms in the potential equation are approxi- 
mated by their first-order expressions. 

This method transforms the nonlinear potential equa- 
tion into an inhomogeneous, linear problem. 

Iteration procedures have been successfully applied 
for the steady-state case by Lighthill’ and Van Dyke® 
for bodies at an angle of attack. The present analysis 
is an extension of Lighthill’s approach to the unsteady 
case of harmonic pitching. The “hybrid” theory pro- 
posed by Van Dyke® has also been extended to the 
unsteady case for cones by Tobak and Wehrend.’ 
The present approach was first suggested by Miles.® ° 
Comprehensive bibliographies of the steady problem are 
given by Van Dyke® and by Sears and Adams." For 
the unsteady case, the reader should consult references 
2, 3, 7, 8, and 9 for a more comprehensive bibliography. 


(4) Formulation of the Problem 


(4.1) Potential Equation 

It is well known that the second-order flow field 
about a slender body is isentropic and may be derived 
from potential flow.® © The approximations implicit 
in potential flow and the differential equations satisfied 
by the perturbation potential are given by Miles.* ° 
In the present notation we have 


= 2,, + 2V2-V2, + (V@/2)-V(VQ)? (4.1.1) 


where 


a? = a,? — (y — 1){2, + (1/2)[(V2)? — (4.1.2) 
The notation is dimensionless (xv = x'//; y = y'//; 
zs =3'/l; 7 = Ur'/l; a, = 1/M) where / is the body 


, 


length, U is the free-stream velocity, and x’, y’, 2’, 
and +r’ are dimensional space variables and time, 
respectively. We note that the differentiation with 
respect to time is a Newtonian reference frame and Q 
is the total velocity potential—.e., the fluid velocity 


is defined by 
V = UVQ (4.1.3) 
We define the perturbation potential by 
V2 = 1,cosa+ 1l.sna+V® (4.1.4) 


The coordinate system is aligned with the principal 
body axes as shown in Fig. 1. If the angle of attack 
is considered to vary with time, then the partial deriva- 
tive with respect to 7 for small angles of attack as 
measured by an observer moving with the harmonically 
pitching body axis system (x, r, @, ¢) is given by 


[O( )/Or] = ( — ar sin ), + 


a(x — [sin 0( ), + (1/r) cos (4.1.5) 
From Eq. (4.1.4) one may write 
Q = (x — xg) cosa+rsinasinéd+ ® (4.1.6) 


Combining Eqs. (4.1.3)—(4.1.6), the pertinent terms in 
Eq. (4.1.1) become (in terms of body coordinates) 


Q., = Py + a(x — sin 6 + cos + 2a(x — X@) E sin 6 + 
r r 


Po, 
Gr®,, sin @ — a’r(x — sin sin + — cos 


[a(x — sin? + sin cos 6 — {| — } | + —— sin 6 + - D, 
or \r r 00 r f 


(4.1.7) 


vQ2-VQ, = a+ ®,) ( ) + (sin sin 6 + ©,) ()+ 
fa) or 


(sin a cos + »,) + ae xe) | sin + 
r r O86 


and 


bol — 


5(V2-V) (VQ)? = 


cos — ar®, sin (4.1.8) 
r 


(cos a + ®,) ) + (sin sin 6 + &,) ()+ 
x or 


2 
(sin acos + »,) = ( I} Leos a+ ,)? + (sin asin 6 + + (sin a cos 6 + (4.1.9) 


The well-known first-order slender body theory de- 
veloped by Ward‘ for steady flow and Miles! * for un- 
steady flows demonstrates that the first-order slender 
body potential, ¢, satisfies, for small a, 


V22¢ = 0 (4.1.10) 


where V..” is the two-dimensional Laplacian operator 


4.1.11 


vie 
Or? r 


The first-order results are 


ao(x) In r + bo(x) + [ai(x, sin @ (4.1.12a) 


g= 


= g(x, r) + ¢1(x, 7) sin 6 (4.1.12b) 


= RR, (4.1.13a) 
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bo(x) = ao ln — In édé (4.1.13b) 
2 0 


and 


a, = + a(x — (4.1.13¢) 


We represent the second-order potential ® as the sum 
of the first-order solution, ¢, plus a second correction 
potential, ¥. Thus, 


(4.1.14) 


By approximating ® by ¢ in the right-hand side of 
Eq. (4.1.1), utilizing Eq. (4.1.14) for the left-hand side 
of Eq. (4.1.1), and observing Eq. (4.1.10), it can be 
shown that the second-order correction potential, VY, 
must satisfy the following equation for small a: 


= By, + M*{2¢n + + 
M 2! 26 sin 0[¢, + (x — XG) X 
(Ger + ¢r¢rr)] + &(x — Xe)¢, Sin 6 + 
2[a sin O( ger + + Frere + + 
+ O(M?a*e?, In? In (4.1.15) 


By using Eqs. (4.1.7) to (4.1.9) and (4.1.12) to (4.1.14) 
one can demonstrate that the neglected terms in Eq. 
(4.1.15) are of the order 


A(V2¥) = O(M*e!, (4.1.16) 


The error estimates will be refined a posteriori, but 
let us defer the error discussion for the present to pur- 
sue the formulation. With the above approximations, 
the nonlinear terms in Eq. (4.1.15) become known 
functions and the second-order differential equation 
becomes the nonhomogeneous two-dimensional Laplace 
(or Poisson) equation in planes transverse to the body 
longitudinal axis. Using Eq. (4.1.12), Eq. (4.1.15) 
can be separated into axial and cross flow components 
and posed in the form 


= Ao(x, r) (4.1.17a) 
= Ay (x, r) (4.1.17b) 
where 
W = W(x, 7) + Vi(x, 7) sin 6 (4.1.17¢) 
Ag(x 
No = Co(x) + In + (x) + (4.1.18a) 
and 
Ai(x, As(x, ¢ 
+ — (x,t) (4.1.18b) 
r r r 
The coefficients are given by 
Co = B°bo,, (4.1.19a) 
ho = (4.1.19b) 
= 2M (4.1.19¢) 


= — M*a,' 


(4.1.19d) 
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Fic. 1. Geometry of the body. 


a} 2(8? + M?)ao,} + 
+ M*)(x — x@)ao, + 
(48? + 6M2)ao} + 
M2[3R? + 5(x — xe)aod} + 
a{ M2R?-(x — x¢)} (4.1.19) 


As = —2M2[al + Rao,} + 
— xg)ao? + R%ao,(x — xe) + 
2R%ay + ao?} + R%ao-(x — xq)}]  (4.1.19f) 


= 402 + a(x —  (4.1.198) 


Eqs. (4.1.17) govern the second-order potential in 
the neighborhood of the body. In order to satisfy the 
boundary conditions at infinity it is necessary to use a 
different approximation to VY which will be called the 
asymptotic solution. Before considering this problem, 
it is necessary to formulate the boundary conditions 
and examine the implications of second-order theory. 


(4.2) Boundary Conditions 


The two well-known fundamental boundary condi- 
tions are (cf. Milne-Thomson!') 


(V2 — V,)-n = 0 (4.2.1) 
on r = R(x), and the evanescence of the disturbance 
potential 


lim (x, 7, 6, t) > 0 (4.2.2) 


For supersonic flow with attached nose shock,*: ® 
r, 0, t) = 7, 0, t) = 0 (4.2.3) 


By applying the relations (4.1.12) to (4.1.14), the 
condition (4.2.1) may be expressed as 


= R, — [a + a(x — sind (4.2.4) 
and 
WV, = Rr¢or — &RR, (4.2.5) 


The conditions at infinity and the initial conditions 
on x are taken care of by observing that ¢ represents 
an approximation for small 7 to the first-order problem. 
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Similarly, WY merges with its asymptotic solution for 
large values of r. These conditions will be developed 
further in Section (5) when the solutions are discussed. 
Before proceeding further with the boundary value 
problem, we must consider the methods of utilizing 
these results to obtain the lift and moment coefficients. 


(4.3) Lift and Moment Expressions 


The lift and pitching moment coefficients are re- 
ferred to the base area, S(1) = ze*/? and the body 
length /. They are obtained by integrating the lifting 
pressure coefficient around the body. Neglecting the 
moment of the axial force, the fundamental relations 
are 


1 2x 
C=- av C,(x, R, 6) R sin (4.3.1a) 
0 0 


1 1 
Cn (x — xo)dx f C,(x, R, sin 
re Jo 0 
(4.3.1b) 
By expanding the pressure coefficient in a Fourier sine 


series 


Cp = C,,(x,r) + C,,(x, r) sin nd (4.3.2) 


n=1 


it can be shown that Eqs. (4.3.1) reduce to 


1 
f RC,,(x, R)dx (4.3.3a) 
0 
1 1 
(x xo) RC,,(x, R)dx (4.3.3b) 
0 
or Cu (xo = 0) + xoCr (4.3.3¢) 


(4.4) Pressure Coefficient 


The exact isentropic pressure equation (in the pres- 
ent notation) is given by®~* 


2 1 ¥/(y¥—-1) 
C, = 1— M°?F it 4.4.1 


(4.4.2) 


where F = (VQ)? + 20, - 1 


and 


(VQ)? = (cos a + ,)? + (sin a sin 6 + #,)? + 
[sin a cos 6 + (1/r)® |? (4.4.3) 
‘ 
20, 2450) sin () + 
ones 


\(x — xg) cosa+rsinasin@+ (4.4.4) 


We next apply the slender body order of magnitude 
estimates of Section (4.1). Noting from Eqs. (4.1.11) 
to (4.1.13) that 


= Ine + ae) 
and [o"( )/or"} 


O(e—") 


one can show approximately that 


This latter estimate will be refined at the end of Section 
(5.1), but will serve for the present purpose. If the 
terms of the order 


O(Meae, In 


are excluded from the analysis, then the expression for 
F can be reduced to 


F = }2[@, + a, sin 6] + 62 + 
+ a(x Xg)®, sin 6] + 
— 2ar®, sin 6} 


M*éa, In €) (4.4.5) 
We then expand F in a Fourier sine series 
F = r) + r,t) sin@ +... (4.4.6) 


The coefficients /) and /) are obtained by resolving ® 
into its axial and cross flow components [with the help 
of Eqs. (4.1.12b) and (4.1.17c)], substituting the re- 
sults into Eq. (4.4.5) and comparing the coefficients of 
sin 26 with those Eq. (4.4.6). Carrying out these opera- 
tions yields 


ky = 2(¢0, WVo,) + 


(¢o, + Vo,)¢0, + + O(M% In (4.4.7) 
t) = 2} (¢1, + WV;,) + a(¢o, Wo,) 
¢o,Vi, + ¢1,Vor $0;¥1, + 
a[(x Xe) (¢o, + Wo,) + 

+ Vi.) + 2¢0,¢1,} + O(M2a,? In (4.4.8) 


Expanding Eq. (4.4.1) by the binominal theorem 
and utilizing Eq. (4.4.6) one obtains 


M°F, 
=n 


M*Fy\ . 
>) sin 6+... (4.4.9) 


whence (4.4.10) 


By introducing Eqs. (4.4.7) and (4.4.8) into Eq. 
(4.4.10), utilizing the boundary conditions, Eqs. (4.2.4) 
and (4.2.5), and evaluating the pressure coefficient on 
the body, one obtains the second-order lifting pressure 
coefficient, |,-r, 


C,, + AC,, (4.4.11a) 


where Cy = —2 + rar (4.4.11b) 


AC,, = at —2[(%,, + Vi, + 
2Ri¢0, + 2Rz*] + M*[(2¢0, + + 
2Rz)} nar + &| + Vi, + 
2(¢o, + (x — + 
M?[(2¢0, + Rz*)(¢1,, + +(x — Xe) + 
R}} + 0(Mtaeé Ine) (4.4.1 1c) 
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Having developed the pressure and lift expressions in 
terms of the first- and second-order velocity potential, 
we are ready to consider the solution of the second- 
order potential equation. 


(5) Solution of the Second-Order Potential 
Equation 

(5.1) Solution in the Neighborhood of the Body 

To solve the second-order potential equation in the 
neighborhood of the body, it is convenient to express 
the second-order correction potential, VY, as the sum 
of a particular solution, y, and a complementary solu- 
tion x. One has 
(5.1.1a) 


(5.1.1b) 


Wo Wo + Xo 
W=n+x 


Introducing Eq. (5.1.1) into Eq. (4.1.17) yields two 
sets of differential equations for Yo, Yi, xo, and x1, Vviz., 


Yor + (1/r) Po, = Ao (5.1.2) 
(l/r) yn, = Ay (5.1.3) 
and + (1/r)xo, = 0 (5.1.4) 


+ (1/r)x1, — (1/r?)x1 = 0 (5.1.5) 


With the aid of Eq. (5.1.1), the boundary condition 
on the body, Eq. (4.2.5), transforms to 


X0, = $0, Rr Yo, (5.1.6) 
x1, — &RR, — Vi, (5.1.7) 


By the method of variation of parameters one may 
construct the following particular integrals: 


dor? de 
Yo — 1) + (nr)? + 
Cor? 
4 + Do In (5.1.8) 
= 1 As 1 1 
D 
Gis 


which may be verified by substitution in Eqs. (5.1.2) 
and (5.1.3). 

Now x satisfies the two-dimensional Laplacian equa- 
tion V2"x = 0; therefore, by analogy to the first-order 
slender body solution one can write* ® § 


xo = Ao(x) Ir r + Bo(x) 
xu = A,(x, t)/r 
Applying Eqs. (5.1.6) to (5.1.11), the coefficients 
Ay + Doand A; + D; become 
(Ao -f- Do) = apo, In R -— 


cert 
4 (2InR—1)+A2InR oR? 


(5.1.10) 
(5.1.11) 


In x + ap In dy"(x — (5.1.12) 
0 


JOURNAL OF THE AEROSPACE 


SCIENCES—OCTOBER, 1960 


and 


> 


MR? 
(A, + Di) = + (1 + 2In R) 


rs (1 3 
( —inR)-- — 
2\2 8 R 


Substituting Eqs. (5.1.12) and (5.1.13) into Eqs. 
(5.1.8)-(5.1.11) and then putting the results in Eq. 
(5.1.1), one obtains 


— a, + ER? (5.1.13) 


2 i. 
W(x, 7) = (nr — 1) + = (In vr)? + 


as + + (Ao + Do)In r + Bo(x) (5.1.14) 


rinr A3 
Wi(x, t) = = > In r + + 


+ (Ai + + Er (5.1.15) 
8r 


It is seen that the coefficients E and By cannot be 
evaluated by considerations of the boundary condi- 
tions on the body alone. Therefore, it is necessary 
to consider the second-order representation of the 
asymptotic solution so that boundary conditions at 
infinity and the initial conditions on x may be invoked. 
It will be seen that after the asymptotic solution is ob- 
tained, the coefficients E and By are obtained by con- 
sidering the behavior of the solution for small values 
of r. 

Before passing to the asymptotic solution, we note 
a posteriori that more accurate order of magnitude 
estimates for Vo and VY, are 


Wo = In? e) (5.1.16) 


and WV, = In (5.1.17) 


(5.2) Asymptotic Behavior of the Solution 


In order to evaluate the functions By and E one may 
follow the procedure of Lighthill’ and require that the 
second-order solution must be compatible with the 


. linearized solution—i.e., that in which the terms of 


O(¢?) are neglected, but those of O(ag) are retained. 
It will be seen that when the behavior of the asymptotic 
cross flow solution is examined for small 7, then there 
arises a term proportional to 7. The coefficient of 
this term is identified as E. With regard to the axial 
flow solution, it will be seen that the additive function 
of x reduces to d(x), and, hence, the second-order 
correction, Bo, vanishes. The procedure adopted here 
renders the asymptotic solution compatible with the 
solution in the neighborhood of the body, so that the 
asymptotic solution describes the flow at large dis- 
tances from the body, and the solution derived in Sec- 
tion (5.1) describes the near field behavior. The fol- 
lowing is a brief outline of the analysis. 

The asymptotic solution will be designated as ¢ and is 


and 


and 


Th 
and ( 


wher 
in se 


= 
en 
sepi 
defi 
mot 
F 
of | 
and pow 
_ _ _ 
a 
| 
It 
with 
the t 
Ther 
| 


13) 


Eq. 


14) 


is 


SECOND-ORDER THEORY 735 


separated into its axial and cross flow components by 
defining 


= do(x, 7) + 7) sin 0 (32.1) 


where k = w//U is the reduced frequency for harmonic 
motion. 
From Eq. (4.1.15), the linearized approximation is 
V2" = + M*\ 26,1 + due} +> 
sin 6[¢, + (x — + 
a(x — xg)¢, sin 0 + 2a¢,, sin 6} + 
O(M?a®, M*g*) (5.2.2 
Putting Eq. (5.2.1) in Eq. (5.2.2), neglecting terms 
of O(M°a*, M*g*), and equating coefficients of like 
powers of sin 6 yields 


= Bdo,, (5.2.3) 
= Bgi,, + = + 
M?{2ik[0, + (x — — 
— + 2¢0,,} (5.2.4) 


From Eqs. (4.2.3) and (4.2.4) the first-order axial 
and cross flow boundary conditions are 


$0, = (5.2.5) 
= —[1 + tk(x — x¢)] (5.2.6) 


lim (do, 0 (5.2.7) 


and the initial conditions 


go(0, r) = (0, r) =0 (5.2.8) 
r) = gi,(0, r) = 0 (5.2.9) 


It is expedient to apply the Laplace transformation 
with respect to x on Egs. (5.2.3) and (5.2.4). Denoting 
the transformed variables by bars, one has 


ds, = o(x, r)} = e—*o(x, r)dx (5.2.10) 
0 


Then Eqs. (5.2.3) and (5.2.4) become 
(V2? — Bs*)Go = 0 
(V2? — = 2iksdi — + 


Ss 


(5.2.11) 


E (Go,) + + (5.2.12) 


The solution to Eq. (5.2.11) satisfying Eqs. (5.2.7) 
and (5.2.8) is 


bo = —a(s)Ko(2) 


where zg = @sr. For small z, Ko(z) may be expanded 
in series yielding 


& -a)| - y¥+1n x 


(1 += *) | = 


4 


(7 + In *) sa(s) + @s) nr — 


Bsr 
a(s) E (7 + In (5.2.13) 


By taking the inverse transform of Eq. (5.2.13) and 
comparing the leading terms with the first-order axial 
flow solution, one finds that 


a(s) = Lfao(x)} (5.2.14) 


and 
ji + In sa(s)' = da In 
ls 2 2 


do’ (x — &) In =bo(x) (5.2.15) 
0 


Since the second-order 69(x) is identical to the first- 
order value, then the second-order correction Bo(x) is 
identically zero. We now consider the cross flow solu- 
tion. 

It is convenient to separate ¢; into real and imaginary 
parts by putting 


od: = + (5.2.16) 


If we put Eq. (5.2.16) for ¢; in Eq. (5.2.12), neglect 
terms of order k? (for the purpose of studying thickness 
effects at low reduced frequencies), and put r = 2/8s 
we obtain 


B°s°V = s) + 2M*sdo, 
= B*s*5(s, 2) + 
sti + bo, — (d/ds)(sdo,) — sxedo,} (5.2.18) 
2 10 1 


where v.22 = — —-—— — 
Oz? z Oz 2? 


(5.2.17) 


(5.2.19) 


In solving Eq. (5.2.18) we replace a by its first-order 
or homogeneous solution, /, which satisfies 


V.7f = (5.2.20a) 
and the boundary conditions 
lim f(s, r) > 0 (5.2.20b) 
lim f(s, r) ~ r)]} (5.2.20c) 
r—0* 


The last equation forces agreement with the slender 
body solution for small r. 
The solution to the problem posed by Eqs. (5.2.19) 
and (5.2.20) is 
f(s, r) = Ki(Bsr) (5.2. 
where = B&,{ 2ao(x)} 


Then using Eq. (5.2.21), Eq. (5.2.18) becomes an in- 
homogeneous problem whose solution can be repre- 
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sented as the sum of a particular solution vp plus a 
homogeneous, complementary solution vy = g. The 
particular solution represents the asymptotic second- 
order solution, while the homogeneous solution repre- 
sents the asymptotic first-order solution. As in the 
case of the real part, we force the homogeneous solution 
to agree with the imaginary part of the first-order 
slender body solution for small r; viz., 


lim (kg(x, 7) > Im| ¢1, r)} (5.2.23) 


r—>0* 


We then pose the solution to the homogeneous 
problem as 
g(s, r) = B,(s)Ki(z) (5.2.24) 
2.23), we obtain 


2ao(x)(x — 


where, by invoking Eq. (5. 
B,(s) = B&,{R? + (5.2.25) 


and we have again used the evanescence condition, 
2(s,r) = Oasr— We consider now the particular 
solution. Substituting = —d(s)Ao(z) and Eq. 
(5.2.21) for uy = f into Eqs. (5.2.17) and (5.2.18) yields 


(V22 — 1)a = (2M?2/8)a(s)Ki(2) 


X 
Bsa(s)zKo(z) + Pi(s)Ki(z)} 


(5.2.26) 


(V2? — = 
(5.2.27) 


where 
P,(s) = {Bs[2a(s) — sxga@(s) — (sd/ds)a(s)]} (5.2.28) 


The above equations represent the nonhomogeneous 
problem to be solved. Eq. (5.2.26) is identical to the 
one solved by Lighthill for the steady-state case. The 
particular solution to @ satisfying the conditions at 
infinity is 


—2M*a(s 
ii(s, z) = |= 3 @) [ + 
his) (5.2.29) 
By the method of variation of parameters, one can 
show that the solution to é satisfying the conditions at 
infinity is 
2M? { 
f Li(m) Ki(n)ndn + 


hie) (5.2.30) 


The solutions (5.2.29) and (5.2.30) are ‘‘exact”’ (with 
regard to the low- Siemans asymptotic problem). The 
problem now is to expand the solutions in series for small 
values of r and determine the coefficient F, as discussed 
in Section (5.1). For the purposes of Section (5.3) we 


1960 


“OCTOBER, 


shall define 


(5.2.31) 
where ¢ and j are the first and second bracketed 
expressions, respectively, in (5.2.30). Finally, we note 


that @ = 


(5.3) Determination of E(x) From the Asymptotic 
Solution 


To begin the discussion, it is necessary to observe 
that in the neighborhood of the body, the second-order 
potential is, from Eq. (4.1.17), 


WV = W(x, 7) + Vilx, r) sin 6 


where WV possesses a term, r(x). It has been already 
noted that / cannot be determined from boundary 
conditions on the body. The second-order asymptotic 
solution has been defined by Eq. (5.2.1) as 


= d(x, r) + 7) sin 

whence, by Eq. (5.2.16), 
= go + alu(x, r) + tkv(x, r)] sin 6 | 
= gd + [au + a] sin 6 


Now, for small 7, « and v have terms of the form 
E,(x)r and E,(x)r where 


(5.3.1) 


E, = (0/0a)E(x, a, &) (5.3.2) 
and E, = (0/0a)E(x, a, &) (5.3.3) 


The solution to @ has already been identified with the 
steady-state solution of Lighthill.6 The small argu- 
ment approximation given by Lighthill is (as applied 
to 7) 


_ MP 2 
im f(z (in (5.3.4) 


lim f(z) ~ 


The procedure for evaluating ¢ is identical to that 
utilized by Lighthill. The analysis consists of inte- 
grating the second term by parts and utilizing the small 
argument representations of the modified Bessel func- 
tions to handle the first integral. 

The integration by parts of the second integral in ¢ 
can be accomplished with the aid of the relation® 

f = 1/2 (5.3.5) 

The details of the analysis are given by the author.'” 

The end result is 


lim ¢(z) ~ 
Bs 


Placing z = @sr in Eqs. (5.3.4) and (5.3.6), factoring 
out the coefficients of 7, one obtains 


M?P, 
E,(s) = {as E (in +y+I1n s)}} 
(5.3.7) 


+ 0(z* In? (5.3.6) 


and Ex(s) = 


(5.3.8) 
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The contribution of vy = g comes from putting z = 
3sr in Eq. (5.2.24) and utilizing Eq. (5.2.25). The re- 
sult is 


Evy(s) = 
Prt 1 
B,(s) E (7 + In : + In (5.3.9) 


Taking the inverse Laplace transforms of (5.3.7)- 
(5.3.9) and adding up the terms then, in the x domain, 
E(x) becomes 

M?ay(x) 


4 


Es) = 


B° (1m ; ) [2ao(x) + (x — xg@)ao'(x)] — 


2 
0 


(x — X@) ao” (x — In (5.3.10) 
0 $ 


From Eq. (5.2.28), the inverse transform of P;(s) 
becomes 


2Bsa(s) — Bs*xga(s) — 
Bs*(d/ds)a(s) } 
(5.3.11) 


II 


B} 2ay'(x) — x@do"(x) + | 
(d?/dx*) [xao(x)] } 


Eq. (5.3.10) represents the term needed for the un- 
steady case [by virtue of Eq. (5.3.3)|. The inclusion 
of this term makes the second-order slender body solu- 
tion in the neighborhood of the body compatible with 
the boundary conditions at infinity. Eq. (5.1.15) 
already includes the effect of the coefficient / on the 
boundary conditions on the body; therefore, Light- 
hill’s method has been extended to solve the problem 
of oscillating body of revolution for low reduced fre- 
quencies. The next section considers the application 
of the second-order corrections to unsteady lift and 
moment coefficients. 


(6) Computation of Lift and Moment 
Coefficients 


(6.1) Limiting Cases 

The unsteady contribution to lift and moment can be 
computed from Eq. (4.3.1) using the second-order pres- 
sure coefficient expression (4.4.11). The pertinent 
velocity potential expressions are Eqs. (4.1.12) and 
(4.1.13) for the ¢ terms and Eq. (5.1.15) for ¥; utilizing 
Eq. (5.3.10) for the evaluation of Ea = E,. It can be 
shown that the present analysis gives results identical 
to those of Lighthill for the steady lift and moment for 
a cone—for example,° 


3 
= — + — — In — 6.1.1 
Cm, la=0 3 haw 6 (6.1.1) 


The first term represents the first-order slender body 
result. The above result can be obtained by following 
the calculation scheme outlined above and utilizing the 
body profile equation for a cone—i.e., 


R=e (6.1.2) 
One also needs the coefficient, L, = E, of the pre- 
vious section. Lighthill shows that this term is given 
by® 
(x 
B, = {- + 


(2M? — 1) ist (6.1.3) 
2] o (: n 6.1. 


No further details of the computation for the steady 
case will be discussed. Now consider the unsteady 
case. 

It can be shown that the inclusion of the first two 
terms of Eq. (4.4.11); 1.e., 


9) 
Cy? = —2 + rar 


leads to the well-known first-order slender body results 
(cf. Miles,': § Dorrance,’? or Sacks*). For the harmonic 
pitching case one obtains 


(6.1.4) 


Cr, = 2[1 + — (6.1.5) 
Cr, = + (1 — (6.1.6) 
and = —2{(1 — xe — Qo) + 
— 2xeQ1 + x6?Qul} (6.1.7) 
Ca = —2(1 — (6.1.8) 


where the moments are referred to the axis of rotation, 


n S(x) 
2, = x dx 


where .S(x) is the local body cross-sectional area and S(1) 
is the base area (all dimensionless). 

The remainder of Eq. (4.4.11) represents the second- 
order correction to the lifting pressure, and the lift and 
moment corrections are derived by carrying out the 
indicated integrations. 

The solution for constant vertical acceleration, 
(Cr, and Cug,.,) is obtained by noting that 


and 


(6.1.9) 


(6.1.10) 


lim [—k’xgay| = w/U = ava 
where a, refers to the solution for harmonic variation 
ofa. Then 

ra) fa) 
——(C, ,Ch = — Cy 6.1.11 
(tay? Cray) = (Co En) (61-11) 
For the first-order case, Eqs. (6.1.5) and (6.1.6) yield 
the well-known results? 7 


Cr. = — 2Qo (6.1.12) 
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(6.2) Second-Order Effects 


TABLE 1 
The present analysis provides a tool with which the Computation Sequence 
second-order thickness and Mach number effects on Terms Equations 
dynamic stability derivatives can be determined for a g and derivatives (4.1.12) 
ge ado, bo, ay (4.1. 13a, b, c) 

body of arbitrary meridian profile shape. For a body VY, (5.1.15) 
whose profile can be described by a polynominal, the As As f, g) 
integrations are quite straightforward, though lengthy. Pix) (5.3.11) 

As an example, the simplest case will be chosen, that of AC, (4.4. 11c) 

ACz, ACn (4.3.3a, b) 


acone. The cone case also provides a comparison with 
the work of Tobak and Wehrend’ who have extended 
Van Dyke’s “‘hybrid theory’’® to the unsteady case. 
This comparison will be discussed later. 

If one defines the difference between the first- and 
second-order lifting pressure as AC,, and separates the 
contributions of the E term from the remainder of Eq. 
(4.4.11), one obtains, for the unsteady lifting pressure 


AC), = AC,,, + AC (6.2.1) 


Pipe 
Table 1 demonstrates the sequence of substitutions. 
In carrying out the above scheme one needs the partial 
derivatives of ¥; with respect to x and /¢, evaluated on 

= R(x) as well as the derivatives of E, with respect 
to x. The general formulas are given below in the 


msanteiatin form of the partial derivatives of these quantities with 
respect to a and a. 
3 
+ — + &(Rao,, + 2RAo,) + 2(RE, + ER;) (6.2.2a, b) 
OR 
here 


af + M*)ao, — k2M?[3R? + 5a-(x — + + + + M*)ao,-(x — — 
k?M?R?-(x — (6.2.3a) 
= af —2M?[3ao? + Rao, — — + —2M?[2aoR? + ao? + (x — xe) + 
— xg)]} (6.2.3b) 
As = af + — (6.2.3¢) 
a{2(B? + M*)ao,, — + 5ao,-(x — + 
(48? + 6M2)ao, + 2(8? + + do,,-(x — — + 2ao-(x — x¢)]} (6.2.3) 
= al —2M?[Sacao, + R%ao,, — k?{aoR? + 2a0?-(x — xg) + 
— xa) + —2M2[6ao? + 3R%ao, + 2acao, + (x — Xe) + 
do,,R?+(x — xq)]} (6.2.3) 


2 
II 


I 


= af 4M? + + — xg) + 2R%ado,-(x — xg) + (6.2.3f) 
+ 41,00,) = af + &{ + ao,R? + 4aodo,-(x — (6.2.3g) 

(Ws, Vs.) = (2. =) E (1 + 2InR)M, — oR ~ + &Ra, — 2RE | (6.2.4a, b) 


= &{2(B? + M*)ao, — + 5ao-(x — x@)]} — 
Raf [(48? + 6M2)ay + 2(8? + M2)ao,-(x — — k®[M2R?-(x — xe)]}  (6.2.5a) 


As, = —2M4{ + R%ao,] — k?[aoR?-(x — xe)]} + 
2aoR? + ao? + 2ao?+(x — xg) + R%ao,-(x — xX@)}  (6.2.5b) 


= — k?a(x — (6.2.5c) 
= — k?a[2R,ao-(x — xg) + (6.2.5d) 


The contribution from the E term can be expressed 
as 


(4.3.3a) and (4.3.3b) and the pressure coefficients 
computed by the scheme outlined in Table 1. 
Numerical values for CL, = and Cu, Mo 


dAC,, /da = —4{E,R + RE,, + E,R,} (6.2.6) 


Pip 


The stability derivatives are computed using Eqs. 
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have been computed and plotted in Figs. 2 and 3 for 
the case of a cone at low reduced frequency, with the 
pitching axis, xg, and the reference axis, %, both taken 
as zero. It is to be noted that the term Cus, for the 


harmonic pitching case as derived from the present 
theory is equivalent to Cu, “a where dy, is the con- 
stant vertical acceleration, and g the constant pitching 
velocity according to conventional stability notation. 
(Note that all of the above terms are in terms of the 
dimensionless variables.) 

One observes that the effects of increasing thickness 
and Mach number are both destabilizing as regards to 
damping moment and cause a substantial reduction 
in the unsteady lift force as compared to first-order 
slender body theory (« = 0° line). 

The effect of pitching axis location on damping 
moment and C,, for a 10° cone is shown in Figs. 4 and 
5. Calculations show that the predominant second- 
order contribution to OCn, ,/OXG (i.e., due to the axis 
shift) is from the reduction in lift, although there is a 
small contribution due to the change in the lift with 
axis location. 

We observe from Fig. 4 that one always has a finite 
damping moment for pitching axis location through the 
base according to second-order theory, and that there 
is actually a crossover between first- and second-order 
theory between 0.75x,g and 0.80xg. We find for 7 < 
1.5 that at the cone centroid axis, x, = 0.75, the first- 
and second-order damping moments are about equal. 
This, however, merely represents the crossover of two 
divergent curves. As the axis moves forward the 
second-order damping moments are of the order of half 
those predicted by first-order slender body theory. 


(7) Discussion of Results 


It is seen that for low supersonic Mach numbers, the 
second-order thickness effects produce a reduction in 
both the out of phase lift and the damping moment 
coefficient (Cur, < 0 is positive damping). For 


a given thickness ratio, the lift and moment coefficients 
both decrease with increasing Mach number, and, for 
a given Mach number, the coefficients decrease with 
increasing thickness. The above trends are true except 
for values of the hypersonic thickness similarity 
parameter, A; greater than about 0.7, where’* 


K; = M6 (7.1) 


This divergence of the trends is similar to the behavior 
of Lighthill’s steady-state solution when compared with 
more exact steady-state theories of Van Dyke® and the 
numerical studies based on the method of character- 
istics presented in the M.I.T. tables for yawed cones.'4 
The point of reversal of the trends is probably the limit 
of applicability of the present second-order slender 
body method. 

Recently, Tobak and Wehrend’ have extended the 
“hybrid”’ theory of Van Dyke’ to the case of unsteady 
flow for constant vertical acceleration, ay,4, and con- 


25 


SHock DeTACHMENT 


10 45 20 25 3.0 
MACH NUMBER 


Fic. 2. Second-order unsteady lift coefficient due to harmonic 
pitching about the axis x; = x9 = 0 at low reduced frequency 
vs. Mach number for cones. 


stant pitching velocity, g. This method combines the 
first-order cross flow solution with the second-order 
axial solution previously derived by Van Dyke,® as well 
as using the exact isentropic pressure equation. The 
sum of the pitching and vertical acceleration solutions 
of Tobak and Wehrend are equivalent to the harmonic 
& contributions of the present theory. Figs. 2 and 3 
provide a comparison of the present results and those 
of Tobak and Wehrend for the case of a 10° semiapex 
angle cone. The present theory is seen to give some- 
what larger reductions in the coefficients for a given 
K, than does reference 7. 

It is anticipated that the role of the present analysis 
will be to provide an explicit set of formulas for com- 


puting C,, ,, and Cy, ,, for bodies of arbitrary 
VA VA 


profile shape for 1/7; < 0.7. 
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coefficient for a 10° semi-apex angle cone, oscillating harmonically 
in pitch at low reduced frequency. 


Some estimate of reduced frequency effects up to 
k® < 1 could be obtained, although the asymptotic 
solution would be approximate. 

Finally, we remark that second-order solution due 
to vertical acceleration could be obtained from the 
present harmonic case by the procedure described by 
Eqs. (6.1.10) and (6.1.11). 


(8) Concluding Remarks 


The present work provides an extension to the un- 
steady case of the second-order slender body theory for 
bodies of revolution developed by Lighthill.6 The 
analysis, in principle, provides an estimate of the effect 
on aerodynamic stability derivatives, of body thickness 
ratio, Mach number, pitching axis location, and re- 
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duced frequency for bodies of revolution with arbitrary 
meridian profile. 

Specific numerical results have been computed for the 
case of a cone with pitching axis at the nose at low re- 
duced frequency. These results indicate reduction in 
the unsteady lift and pitching moment. The reduc- 
tions increase with both Mach number and body thick- 
ness ratio. There is, however, an apparent reversal 
of trends beyond a hypersonic thickness similarity 
parameter, J/;, of about 0.7. This trend reversal 
is considered by the author to indicate a breakdown of 
the theory rather than having any real significance. 
Comparison of Lighthill’s steady-state results’ with 
more exact results (cf. Van Dyke®) seem to bear out 
this conjecture. 

It is believed that the most important conclusion of 
the present paper is the demonstration of the impor- 
tance of the second-order cross flow correction in the 
unsteady lifting problem for low supersonic Mach 
numbers. 

It is possible that some of the limitations in the present 
theory could be removed by the use of a more exact 
first-order theory in the iteration scheme. In_ this 
connection the work of Dorrance? might provide a 
starting point. 
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On the Response of the Laminar Boundary 
Laver to Small Fluctuations of the Free-Stream 
Velocity! 


NICHOLAS ROTT* ann MARTIN L. ROSENZWEIG** 
Cornell University 


Summary 


The linearized treatment of small time-dependent disturbances 
of a laminar boundary layer, initiated by Lighthill, is extended in 
several ways. In particular, the high-frequency expansion is 
continued beyond the leading (Stokes) term. Several interesting 
questions of ‘“‘joining’’ occur, which are discussed but left unre- 
solved. In addition, a practical method for obtaining the re- 
sponse to the laminar boundary layer to an impulsive change in 
velocity is presented. The methods are applied to the case in 
which the basic steady flow belongs to the Falkner and Skan 
family of similarity solutions. 


(1) Introduction 


IME-DEPENDENT boundary-layer flow problems 
gp pee in connection with many interesting and im- 
portant fluid mechanical effects. Many of these phe- 
nomena, however, are not even qualitatively explained 
by present theories, and a long-range effort in the 
theory will be needed to bridge this gap. 

The problem of the response of an otherwise steady 
laminar boundary layer to small disturbances has been 
recently treated in several papers. It is intended here 
to present a survey of this linearized theory and to 
generalize and extend it in several ways. A practical 
method of computation for arbitrary free-stream fluc- 
tuations is one of the results. The value of these inves- 
tigations in closing the forementioned gap is still in- 
direct. The fundamental question in most cases in- 
volving time-dependent viscous phenomena is: to what 
extent can a linear theory describe effects which are 
basically nonlinear? Even to answer this question it 
is clear that a well founded, rather complete and prac- 
tical linear theory is a necessary tool. 

In Section (2) of this paper, the general linearized 
problem, as first set forth by Lighthill, is described. 
Succeeding sections treat the low- and high-frequency 
solutions, especially as applied to the Falkner and Skan 
family of similarity solutions. In Section (5), the 
response of the laminar boundary layer to a small im- 
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pulsive change in the free-stream velocity is discussed, 
and a practical method of computing this response is 
presented. 


(2) Description of the Linearized Problem 


Consider a free-stream velocity U(x, t) which con- 
sists of a basic steady distribution l(x) with small 
periodic fluctuations of amplitude el\(x) superimposed 


U(x, t) = Uo(w) + (e<1) (1) 


Using Eq. (1) as a free-stream boundary condition, 
Lighthill' treats the “‘linearized’”’ problem for incom- 
pressible, laminar, two-dimensional flow, by assuming 
a solution to the boundary-layer equations of the form 


~ 


u = U(x, ¥) + (x, 
v = w(x, ¥) + 


a 


where m% and vw are the velocity components in the 
steady case. Retaining terms of order ¢ only, ™ and 2, 
fulfill the equations 


Ou; Ov 
Ox oy 
4 Ou; 4 Ou 4 Ou; Ou 
) o— — — 
: 
oy? 
u(0) = (0) = 0, m(o) = U, 


Lighthill proceeds to solve these equations by making 
two separate approaches, for low and high values of the 
frequency w, respectively. 

For small frequencies the series expansion 


@ @ 
u = >> 11 = (iw) (4) 


n=0 n=0 


is introduced. For the special case when the amplitude 
of the fluctuating free-stream velocity is proportional 
to the steady value, i.e., for 


U(x, t) = + ee) (la) 


(U, = Up), Lighthill gives the explicit solutions for 
uy and vo, expressed in terms of the steady state 
solutions u% and v% as follows: 
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(5) 


= 2 Vo ay 

These expressions fulfill Eq. (3) for w = 0 and can be ob- 
tained by differentiating m and v with respect to the 
free-stream velocity U(x), a procedure which is easily 
understood by the interpretation that 1» and v4 are due 
to a ‘‘quasi-steady”’ change of to Uo(x)(1 + 
Lighthill bases the determination of “, (for the special 
case U,; = Uy) on the Karman-Pohlhausen method and 
proposes after many further simplifications (justified 
by rather lengthy error estimates) the expression 


Ub? 
= - (i Ane + (2.1 (6) 


where 6 is the Pohlhausen boundary-layer thickness, 
np = y/6, and with the Pohlhausen parameter A = 


A = (36 — A)/240 (7) 
The shape of the profile (6) is chosen such that 
Uo 


at y = O for all values of A. The further significant 
properties of the profile given by Eq. (6) are the fol- 
lowing: 

(i) If both the velocities 1, and the frequency w in 
Eq. (4) are to be made dimensionless, Eq. (6) indicates 
that the proper ‘‘frequency parameter’”’ k (say) is 


k = w6?/p (8) 


(ii) The shape of the m, profile depends only on the 
shape of the local steady profile, and its scale is given 
by the local steady boundary-layer thickness. 

As the order of magnitude of 6 in steady flow is given 
by 


5° ~ vx/ Uy 
an alternate possibility for Eq. (8) is to use 
k’ = wx/Up (9) 


For similar steady (Falkner-Skan) profiles the two 
definitions agree for all x except for a constant factor. 
It will be shown that if the fundamental solution belongs 
to the Falkner-Skan family, the expansion (4) leads to 
similar solutions 1, if the frequency parameter given by 
Eq. (9) is used instead of w. This result permits a re- 
interpretation of the approximations made by Lighthill 
to obtain the expression (6) for uw. It can be said that 
4; Will be ‘‘of the form’”’ given by Eq. (6) if the actual 
basic steady boundary-layer development is locally 
“fitted” by a similar steady flow. In the light of 
Lighthill’s error estimates this procedure is justified. 

In the next section, the function 1, will be determined 
for a ‘‘similar’’ basic flow, in which case 1; will have ex- 
actly the properties (but not exactly the shape) given 
by Eq. (6). The exact shapes will be determined and 
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compared with the approximate profiles given by Eq. 


(6). 


(3) Low-Frequency Disturbances of the Steady 
Similar Solution 


It is seen that the use of the frequency parameter (9) 
implies an expansion in powers of a quantity which is a 
function of x. (The only exception is the case L,/x = 
const.; i.e., stagnation-point flow.) For steady flows 
of the Falkner-Skan type (i.e., U) = cx”) it is conven- 
ient to introduce the quantity 

2iwx 


0 
+ m) Uo (1 + m)cx™—! ) 


Besides £, the profiles should depend only on the 
Falkner-Skan variable 


m + 1 Uo y 
= 
ve 
The basic flow is well known: 
Uo F’(n) (12) 
% = — )’F + F’ 
where F represents the Hartree functions. Let UU be 
of the form 
U; = Ox" (13) 


where m, can be different from m, and put 


= Unf, 2) \ (14) 
— (U160)'f + U160'nf, U fe 


where ¢’ = dé/dx. Continuity is fulfilled; introduction 
of these expressions in Eq. (3) leads, after same calcula- 
tions, to the equation 


+ Fas = + (8 + Bi) + 
+h OFS — 2 — — + 
= 0 (15) 


V1 


where 
B 2m/(1 + m) 
and By = 2m,/(1 + m) 


(It is no misprint that m and not m, appears in the de- 
nominator!) The complication caused by m, # m ap- 
pears small. Note however that (with the only excep- 
tion for 8 = 1, i.e., for a stagnation-point flow as basic 
flow) Eq. (15) contains partial derivatives with respect 
to &. 

In what follows, the discussion will be restricted 
mostly to the case m,; = m, where the quasi-steady solu- 
tion is known. For the low-frequency expansion, put 


Equating powers of &, it is found that gp fulfills Eq. (15) 
for f with = 0. For 6, = 8, the solution is 


= (1/2)(F + nF’), go’ = F’ + (17) 
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in agreement with Lighthill’s formula, Eq. (5). The 
equation for g; is 


gi’ + — + (3 — 28) = go’ — 1 (18) 


This equation was solved numerically at the Cornell 


Computing Center for the flat plate (@ = 0), for 
stagnation-point flow (6 = 1), and for a case of de- 
celerated flow (8 = —0.16). The results are plotted 


in Fig. 1. For comparison, the profiles according to 
Lighthill's approximate formula (6) are also plotted, 
with 6 and A (and 6) computed by the Pohlhausen 
method for the corresponding Falkner-Skan cases. The 
agreement with the exact solution is particularly good 
for 8 = 0; it is less satisfactory for both accelerated and 
decelerated flow. Still, considering the approximate 
nature of Lighthill’s analysis, the range of applicability 
of Eq. (6) is surprisingly good. For the shear, which is 
proportional to g,’’(0), the comparison is found in the 
following table: 


Exact Approximate 
B= 1 0.2854 0.320 
B= 0 0.6000 0.619 
B = —0.16 1.0974 0.875 


For n > 1, g, fulfills the equation (for 8; = 8) 


g, + Fg, — + n(1 — B)|F’g, + 
[1 + 2m(1 — B)) F’'gn = (19) 


No solutions have been calculated so far form > 1. It 
is, however, interesting to note that g, = F’ solves the 
homogeneous part of Eq. (19) for all m (including 0 and 
1), for all 8, and even for all @;. This result will cer- 
tainly be needed for the more general solutions. 

The fact that f fulfills a partial differential equation in 
nand £ does not formally cause any particular complica- 
tion when the expansion (16) is made. However, the 
definition (10) shows that the expansion diverges for 
x > 0, if m > 1, and for x ~ o~, if m < 1; uniform 
validity for all x occurs only for the stagnation-point 
flow, m = 1. It should also be noted that the “‘thick- 
ness’’ of the low-frequency solution is determined by 
the ‘‘basic’’ thickness 6). This fact follows from the 
success of the expansion (16), which uses the Falkner- 
Skan variable (11) based on the fundamental steady 
flow. 

A particular case of interest is obtained by putting 
8 = 1, i.e., by choosing a basic stagnation-point flow, 
and putting 6, = 0, which means that the time-depend- 
ent part has a (spatially) constant free-stream velocity 
amplitude. The superposition of such a time-depend- 
ent flow upon the basic pattern results in a stagnation 
“point”’ which oscillates while the ‘‘intensity’’ (that is, 
the velocity gradient) stays the same. In this case, f 
fulfills the equation 


fann + Fon — (F' + Of, +E+1=0 (20) 


It has been shown by Glauert? that the solution has the 
form 


0 0050 ai g' 02 0 025 05 gf 


Fic. 1. The time-dependent boundary-layer profile com- 
ponent, proportional to the frequency, for the low-frequency 
approximation. Plotted are (solid line) the exact solutions of 
Eg. (18) with the stagnation-point flow (8 = 1), the flat plate 
(8 = 0), and a case of retarded flow (8 = —0.16) given as the 
“‘basic’’ steady boundary layer. The broken lines are the 
approximate profiles calculated from Eq. (6). 


(21) 
where ¢(7) fulfills the homogeneous equation 
¢’’ + Fo’ — (F’ + §)¢ = 0 (22) 
with the boundary conditions 
= 1, =0 (22a) 


The profile ¢ can be interpreted directly as the flow with 
a stagnation-point pattern at rest, where the wall per- 
forms oscillations in its own plane. It is connected 
with the f, profile by simple “‘relativity’’ arguments— 
i.e., change of coordinates. The operation requires 
some care; its detailed presentation can be found in 
Glauert’s work. 

The special case under consideration is interesting as 
it can be shown to be the only member of the family 
described by Eq. (15) for which the quadratic terms in 
€ are identically zero, so that Eqs. (20) and (22) hold 
exactly for any «. In addition, it can be easily verified 
that the solutions of Eqs. (20) and (22) even satisfy the 
full Navier-Stokes equations, just as the basic steady 
stagnation-point flow does. The low-frequency solu- 
tions for Eq. (22) were discussed by Rott* and Glauert ;* 
they exhibit the characteristic properties of the general 
linearized case discussed above. 


(4) High-Frequency Expansion 


In the limit w > ©, the leading term is obtained by 
retaining in Eq. (3) the terms with the factor w together 
with the derivative of the highest order: 


iwi, — = (23) 
The solution will be called 1; : it is 
= U,(x)(1 — e~» Vier) (24) 


i.e., Stokes’ solution is found, with the amplitude de- 
termined by the ‘“‘local” value. 
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Lighthill restricts his discussion for large w to the 
solution (24). For a further expansion, the procedure 
employed by Glauert? and by Rott’ [for the simple 
special case, Eq. (22)| suggests the following series: 


= >> (iw)—"/ (x, 


n=0 


(25) 


n=0 


Y = yViw/y (26) 


ll 


where 


A problem arises at the introduction of these expres- 
sions in Eq. (3), as the steady solution m and z% is given 
in function of the variables x and y, while the pro- 
posed expansion calls for the use of x and Y. The dif- 
ficulty is overcome by a procedure used by Glauert:* 
up and v are expanded in function of y near the wall 
y = 0, and the terms in equal powers of (iw)~'/? are 
gathered after the substitution (26) is made in the wy and 
Up Series. 

It is a fundamental property of the high-frequency 
solution {exemplified in the leading term (24)| that in 
the limit w > © the layer affected by fluctuations be- 
comes arbitrarily small. This serves as a justification 
of Glauert’s procedure, since for w > © the unsteady 
boundary layer is restricted to a small fraction of the 
steady layer near the wall. 

Upon application of the above method it is found 
that the “interaction terms’? involving and in 
Eq. (3) will influence the expansion (25) only from the 
term 2," on; that m"" is zero: and that the solution 
up to and including the term 2° is obtained by solving 
the equation 


— v(O7m/Oy?) = + (UoUi)’ (27) 


The solution which remains finite at y > © and ful- 
fills = 0 is 


= [U, + — (28) 


A second difficulty is encountered now: the boundary 
condition = U is not fulfilled. [For = 
0, the trouble is only delayed to the next term, as will be 
seen presently.]| Glauert treated the special case of 
the homogeneous equation (22), and no such difficulty 
occurred. The first analysis in which a high-frequency 
expansion was made without fulfilling the boundary 
condition at infinity is due to Illingworth,‘ who seems 
not to be concerned with any possible complication 
arising from this fact. If it is recalled that the method 
involves an expansion for small y, it appears to be almost 
accidental that the boundary condition for Y > © can 
be fulfilled term-by-term. 

There is no doubt that for physically meaningful 
solutions the terms containing e” must be rejected; 
this condition supersedes and in general disagrees with 
the “‘old’”’ boundary condition = U;. Further 
investigation will show, however, that the appearance 
of terms which are proportional to Y, Y?, etc., in m 
cannot be avoided. The general term ~ ") will con- 
tain a part which is multiplied by e~ ”, and an additional 
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polynomialin Y. The validity of such a solution needs 
further justification. 

The part multiplied by e~ 
noticeable only within a thickness 6; ~ Vy, w; com- 
pared to the basic flow thickness of the order 6, ~ 
V vx/Up, it has the ratio 


” corresponds to a flow 


61/59 = Uo/wx (29) 


(independent of v!) which becomes small for w > ~, 
The polynomial part which is not multiplied by e~ * 
“protrudes” into the outer layer. Now it will be re- 
called that Y is already infinitely large at the edge of 
the outer layer, y = 6); the order of Y there is the 
reciprocal of the ratio (29). Therefore the expansion 
in ascending powers of Y applies to the inner part of 
the outer layer; nothing is known about the conver- 
gence, but the nondecaying part of the solution can be 
interpreted as the high-frequency influence on the 
outer layer, obtained as an expansion starting from the 
‘interface’ between the two layers at about 6) — 6;, 
which is naturally still near the wall (actually, at VY > 
Qasw—> 

In summary, the existence of a high-frequency 
boundary-layer flow affecting the whole basic boundary- 
layer domain is indicated. The method under consid- 
eration gives only the expansion near the wall, and it is 
remarkable that nowhere is a “degree of freedom” 
found to make use of the boundary condition at infin- 
ity. It is therefore a nontrivial assertion that the ex- 
pansion represents a flow which fulfills the free-stream 
boundary condition. 

In a special case, the assertion can be proved. While 
solving the homogeneous equation (22), Glauert? and 
Rott* found no parts of ¢ which would not vanish ex- 
ponentially for Y —~ o. If, however, Glauert had 
tackled the inhomogeneous equation (20) directly by 
the expansion method, he would have found nonvan- 
ishing parts corresponding to the expansion of F’(n) 
in the solution (21). 

It is conjectured that in all cases, the outer expansion 
belongs to a solution which joins the free-stream condi- 
tion at infinity. It would be sufficient to show that 
such an outer solution exists; however, physical 
arguments of existence are not admissible as long as it 
is not known to what extent the simplified equations 
under consideration represent the physical situation. 

We now proceed to the calculation of the term ™‘°. 
It is found that it suffices to approximate m and uv by 
one term only: 


uo = B(x)y, = —B’(x)y?/2 (30) 


An “interaction term’”’ is found, formed by all terms of 
Eq. (3) which are not contained in Eq. (27), to give the 
inhomogeneous part of the equation for ‘*. For this 
purpose, it is sufficient to consider the interaction ef- 
fect caused by u,, given by Eq. (24), and 7; (to be 
determined by continuity). After some calculation, 
the equation for 1 is found to be 
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oy? 


= Vv + 1) + 


1 
Vv UB" 
The solution is 


l 
-1]+ 


Thus, 7°” contributes a constant term and a term 
proportional to Y to the “‘outer’’ solution. 
For the discussion, it is instructive to introduce 


(33) 
where a(x) and 6(x) are dimensionless and of order 1. 


It is seen that » drops out of the amplitude in (382) if 
This effect can be 


B(x) = a(x)( Uo do) = vx 
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= (iwpp)'/?(U, + 


(iw) Ui’ + (5/8) (84) 


where 7(x) is the shear for steady flow. 

This result is already contained in a thesis prepared 
by Gibson® (under the direction of Prof. C. C. Lin). 
Gibson made contributions to several questions con- 
cerning the time-dependent boundary layer, but the 
most significant results of his thesis are in connection 
with the high-frequency approximation. The exhaus- 
tive treatment of this question is not restricted to the 
linearized theory. Gibson finds the ‘‘acoustic stream- 
ing” (i.e., the change in the average velocity distribu- 
tion due to nonlinear interaction) in the “inner” and 
“outer” part of the boundary layer. 

It is straightforward to apply the previous general 
results to the special case of the Falkner-Skan flow. 
If Eq. (15) is to be solved directly, the appropriate ex- 
pression for / for large & is 


the expression (33) is introduced. seal 

traced back to the independence of the ratio (29) from i (*) SF eo"? £65 (35) 

the value of v. 1w n=0 

The results so far can be summed up in the following ie Chins 

formula for the shear at the wall. Using the high- 

frequency expansion terms up to °°’, the shear 7 = (36) 

OV) 1S n=0 

- The first terms of the series are 

hs = F"(0){—8, + (1 — + Be — (1/8)(3 — 68 — 68,) — (37) 
(1/8)(3 — 68 — 26,)V2e~¥ — (1/12)(1 — 28)V%e7' if 


At this point, a second joining problem of the high- 
frequency solution has to be discussed. As the fre- 
quency parameter § ~ x'~” for the Falkner-Skan case, 
the low- and high-frequency parameter solutions will be 
both found in general for one (arbitrary) frequency 
along the wall. For example, for m < 1 the low-fre- 
quency series gives the expansion beginning from the 
leading edge, while the high-frequency series is valid for 
The high-frequency solution is, however, well 
determined: the leading term as well as all the subse- 
quent steps follow without any ‘“‘influence” from the 
low-frequency expansion, although the latter actually 
represents a flow which is upstream of the high-fre- 
quency part. It is therefore a nontrivial assertion that 
both expansions refer to the same flow. As to “‘physi- 
cal arguments” of ‘‘existence,’’ the statement made 
previously can be repeated: it is not known to what 
extent the present linearization covers the physical 
situation; thus, such arguments are not admissible. 
The question whether the two expansions join is open. 

An exception—already mentioned—is the stagna- 
tion-point flow 8 = 1: £ is independent of x, and the 
solution for any frequency applies along the whole wall. 
Eq. (15) becomes an ordinary differential equation with 
£ appearing only as a parameter. It is remarkable that 


large x. 


formally the series solutions do not show any peculiar- 
ity in the special case 8 = 


1, while physically the stag- 


nation-point case is essentially different from the rest 
of the Falkner-Skan family. 

For an arbitrary basic steady flow, the previous dis- 
cussion of the Falkner-Skan case indicates that blunt- 
nosed bodies with an ordinary stagnation point will be- 
have differently from sharp-edged bodies. In the first 
case, the magnitude of the frequency parameter based 
on stagnation-point boundary-layer thickness can be 
taken as a criterion of whether a low- or high-frequency 
expansion should be used for the whole body. For 
sharp edges, there always will be a low-frequency region 
for convex corners and a high-frequency region for con- 
cave corners; thus, matching problems can arise. 


(5) Small Impulsive Change 


The problem of the response of the laminar boundary 
layer to changes of the free-stream velocity with arbi- 
trary (not necessarily periodic) time dependence was 
treated by Moore,® who proposed an expansion (for the 
special case of the flat plate) in a series with the terms 
xU/U?, x°XU/U*... (U = dU/dt) and ascending powers 
of each of these terms for higher approximations. A 
drawback of Moore’s method is that it presupposes 
that the members of the series form a decreasing se- 
quence, which occurs only in certain phases of quasi- 
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steady motion; within these bounds, however, results 
of great generality are obtained. 

The present approach follows the classical pattern for 
linearized problems: if the periodic solution is known, 
then, by Fourier’s theorem, all solutions can be obtained 
in principle by superposition. The fundamental case is 
obtained by determining the response to a small impul- 
sive (step-function) change. All indicial flow problems 
can subsequently be treated by the Duhamel integral 
method. 

The case under consideration differs from the classical 
problem of the boundary-layer development for a body 
starting from rest, as only a small change in a basic 
steady flow will be considered. Also, Stewartson’s non- 
linear problem of the half-infinite flat plate starting 
from rest,’ which was recently solved by Lam and 
Crocco,’ is beyond the realm of this linearized treat- 
ment. 

The response to an impulsive change of the motion 
for the linear problem was recently treated by Rosen- 
zweig.® His first example deals with the case where the 
periodic solution was previously given by Rott* and 
Glauert:? a basic stagnation-point flow, which is 
steady at infinity, with a time-dependent boundary- 
layer motion caused by the wall, y = 0, moving in its 
own plane. The periodic problem involves the solution 
of Eq. (22) quoted above. As the problem is exactly 
linear in this special case, the Fourier-integral method 
will give the exact solution for the case of the impul- 
sively started wall motion. 

Rosenzweig’s investigation is first restricted to the 
shear. Let the results for the periodic flow case, when 
the wall velocity is U, = const. e’, be summarized as 
follows: the shear is 


T» = U Vapp ¢(E) (38) 
where, putting Uj) = ax, ~ is the frequency parameter 
iw/a and ¢ = —®¢’(0) is to be obtained as a solution of 


Eq. (22). Using the values given by Glauert and 
Rott, the high-frequency dependence of ¢ on & is 


= + 0.4622E-! — 7/7... (39) 
for || > 1, while for |¢| <1, 
¢ = 0.8113 + 0.4932 — 0.09478... (40) 


Now, the shear due to the unit jump U,, = const. 1(¢) 
is of the form 


tw = Uy $(at) (41) 


where, as is well known from the theory of Fourier and 
Laplace transformations, the connection between ¢ and 
¢ is given by Bromwich’s integral: 


or, ¢ is the inverse Laplace transform of g(£)/é. 

As ¢ is known only in the form of the two series (39) 
and (40), the determination of % leads to some difficul- 
ties. First, the use of the series (39) is straightforward: 


term-by-term transformation, with the help of elemen- 
tary formulas, leads to the series 


= (wat)~'/? +0.4622(at) — . . 
(43) 


which is valid for small values of at. It is possible, at 
the expense of some more work, to obtain the whole 
profile development for small times, either by the 
transformation method or by direct expansion for 
small at. The leading term is always Rayleigh’s solu- 
tion, corresponding to the Stokes term which leads the 
high-frequency expansion. 

The difficulty arises when the series (40) is used. It 
can be concluded from this series that ¢(~) = 0.8113, 
but the last two terms do not give directly any further 
information on the behavior of ¢ for large at. The 
reason is that = 0 is a regular point of the function ¢, 
and it is known from the theory of Laplace transforms 
that the regular part of an expansion of the transform 
gives no direct information about the original function. 
As the calculation of the two last terms in the series 
(40) involves by far the most labor for the periodic 
solution, it is particularly disappointing that they are 
useless for the transform method. 

What would be needed, in principle, to determine 
the asymptotic behavior of ¢ for large at? One method 
makes use of the ‘‘eigensolutions” of Eq. (22). These 
are the solutions for the homogeneous boundary condi- 
tions ¢(0) = ¢(©) = 0 replacing the boundary condi- 
tions (22a). In words: the wall is at rest; so is the 
stream at infinity. It can be shown that solutions exist 
only for certain discrete eigenvalues £, which are all 
real and negative, i.e., Eq. (22) with homogeneous 
boundary conditions represents a Sturm-Liouville sys- 
tem. It is also physically satisfactory that the eigen- 
solutions, which are disturbances in the absence of any 
(persisting) external forcing, should have a time-factor 
e with & < 0, indicating stability with respect to dis- 
turbances which depend spatially on y only. 

Use can be made of a series of eigenfunctions for the 
determination of the response to a step-function in wall 
velocity. The necessary steps can be described also 
in the language of the Laplace-transform method—..e., 
by considering the evaluation of the integral (42). The 
following theorem is needed: for an eigenvalue &, 
when the problem with homogeneous boundary condi- 
tions has a solution, the ‘‘inhomogeneous’”’ problem, 
with boundary conditions given by Eq. (22a), has no 
solution, or, to be more exact, its solution becomes 
singular. Thus, ¢(£)/é will have singularities at the 
eigenvalues—all (in general simple poles) situated on 
the negative real axis. The singularity nearest the ori- 
gin, i.e., the smallest eigenvalue lel, determines the 
leading term of the asymptotic expansion for large at. 

In summary, it is seen that the asymptotic behavior 
of ¢ can be determined only by a rather big effort, 
which is perhaps not worthwhile. A different approach 
was made by Rosenzweig*® who noticed that a simple 
expression for ¢(é) fits closely and interpolates success- 
fully both the series (39) for large ¢ and (40) for small 
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£; the expression is 


é+ 1.3164 4 0.4622¢ 4 (44) 
Ve + 2.6328 + 3.860)? 

The first term here is chosen so that the first terms in 
the series (39) and (40) are exactly reproduced; the 
second term of the expression (44) corrects for the exact 
value of the second term in each series. The principle 
of this ‘fitting’ is obvious: it is easily extended to 
further terms and adapted to different series. A lucky 
circumstance permits the check of the accuracy of Eq. 


(44) at the point € = 1, where Rott* found an exact 
solution with ¢(1) = 1.2326; Eq. (44) gives ¢(1) = 
1235. 


With ¢ given by Eq. (44), the determination of ¢ is 
an elementary problem; it is found that 


e(at) = 0.8118 erf[(2.6328at)'”?] + 
+ 0.4622at 3-86at (45) 


This expression is valid for all ¢. The exact asympto- 
tic ‘approach exponent”’ is certainly not contained in 
this solution, but this does not exclude a very ‘‘close 
fit,’’ which is assured even if the very small asymptotic 
“approach” terms are not exact. The value 2(@) it- 
self is naturally exact. 

For the corresponding profiles, two limiting cases are 
known exactly: a Rayleigh profile is found for ¢ > 0, 
and the profile for > © is the quasi-steady solution 
(# = 0 in the low-frequency expansion). For inter- 
mediate times, Rosenzweig’ calculates a mixture of the 
two limiting profiles by superimposing these shapes in 
some ratio and at some scale which is to be determined. 
The time-dependent ratio and scale of these profiles is 
found by fitting the assumed solution to as many bound- 
ary conditions as possible, and in particular by fitting 
it to the shear, which is known to a high degree of ac- 
curacy, so that the profile shapes determined by this 
method certainly represent a good approximation. 

Further details are given in Rosenzweig’s thesis. 
Recently, the same case of the stagnation-point flow 
with an impulsively moving wall was solved by Wat- 
son.'° He applied the Laplace transformation method 
to the determination of a profile series valid for small 
times. This expansion, as was pointed out in the pre- 
vious discussion, does not cause any fundamental dif- 
ficulties. Watson gives a large number of terms 
(based on the high-frequency series given by Glauert?) 
so that even the approach to the asymptotic value at 
t{—» ~ can be rather well determined in this particular 
example. The difficult problem of the asymptotic ex- 
pansion is thereby avoided, but Watson realizes that, 
in general, an impractically big number of terms of the 
“small time’’ expansion would be needed to give good 
results for large ¢. He proposes a Karman-Pohlhausen 
type of approximation for the general case, which he ap- 
plies successfully to the example treated by series. 

Clearly, Rosenzweig’s method gives a practical way 
to obtain the response to the impulsive change of the 
free-stream velocity which is easily adaptable to any 
problem where low- and high-frequency expansions of 
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the shear are available. These expansions were dis- 
cussed in full generality in the previous sections of the 
paper; the only limitations for obtaining the step- 
function response are those which were already found 
for the periodic solution. 

A necessary condition for the possible determination 
of a step-function response by the transform method is 
that a smooth joining between low- and high-frequency 
periodic solutions exists. It was pointed out before 
that, for the general Falkner-Skan case, existence is not 
assured. Conditional, however, to existence, the 
Falkner-Skan flow reacts to a small impulsive changes 
of its steady free-stream value = U»(x)1(t)] witha 
flow field of the form 


uy, Uot y 
= f(—,= 46) 
U,; f ( x 


where f can be determined by the method used in the 
previous example. 

Again, there are no matching problems for the stag- 
nation-point flow, m = 1. The third-order equation 
(15) is for 8 = 1 in general more complicated than the 
particular second-order equation (22). But £ is only a 
parameter, and no “‘binding’’ of neighboring solutions 
by derivatives with respect to £ will make the solution 
obtained at a given value of £ questionable. 

As far as the low- and high-frequency solutions are 
known for arbitrary U(x) (the low-frequency end is 
known in Lighthill’s approximation), the response to 
impulsive change can be found in full generality, fol- 
lowing the steps outlined in our previous example. 
Some further details are given in reference 9. 

The solution of more general indicial problems fol- 
lows by classical methods. For instance, if 


U; Uo(x)1()t", and Uo = cx” 


the response is 


n Uot/x 
My x Uot 4 
—— - — (47) 
U, (= ) f f ( x bo 


In summary, it has been shown that the response of 
a laminar boundary layer in plane incompressible flow 
has a practically computable solution for arbitrary small 
changes of the free-stream velocity. However, several 
results are only conditionally valid because of unsolved 
questions of joining (existence).* 
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(Continued on page 787) 


* NOTE ADDED IN PROOF: After the completion of this paper, 
investigations by S. H. Lam and N. Rott resolved the questions 
formulated above with (almost complete) proofs of joining, also 
substantiated by numerical calculations. A Cornell University 
report on these results will appear shortly. 
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Dynamics and Thermodynamics of Re-Entry’ 


W. H. T. LOH* 
Chance Vought Aircraft 


Summary 


This paper presents a summary of the approximate analytical 
solutions on velocity, altitude, range, deceleration, maximum 
deceleration, time rate and maximum time rate of average and 
local stagnation region heat transferred into vehicle, total heat 
transferred into vehicle, and time of flight for boost glide or 
orbital re-entry into planetary atmosphere for various kinds of 
re-entry trajectories. These trajectories include (1) great circle 
flight at small angles of inclination, (2) great circle flight at large 
angles of inclination, (3) minor circle flight at small angles of 
inclination, (4) minor circle flight at large angles of inclination, 
(5) constant equilibrium temperature flight, (6) constant de- 
celeration flight, (7) constant aerodynamic load factor flight, 
(8) constant rate of descent flight, (9) constant angle of inclination 
flight, and (10) constant rate of heat input flight. 


Foreword 


The equations of motion for boost glide or orbital re-entry 
into planetary atmosphere are approximately written into an 
ordinary differential equation which may be used for great circle 
or minor circle flight at small or large angles of inclination. The 
equation includes many terms, some of which represent the 
gravity force, the lift drag ratio, the minor circle angle, and the 
angle of inclination. If these particular terms are disregarded 
or simplified to those used by others such as Allen and Eggers, 
the differential equation and its solutions are also simplified and 
thus yield the same results as those obtained by Allen, Eggers 
and others. 

A number of solutions are obtained. These include the fol- 
lowing: 


(1) Great circle flight at small angles of inclination. 
(2) Great circle flight at large angles of inclination. 
(3) Minor circle flight at small angles of inclination. 
(4) Minor circle flight at large angles of inclination. 
(5) Constant equilibrium temperature flight. 
(6) Constant deceleration flight. 
(7) Constant aerodynamic load flight. 
(8) Constant descent flight. 
(9) Constant angle of inclination flight. 

(10) Constant heat-input flight. 


In most cases, the solution on instantaneous and final velocity, 
density, altitude, distance along flight path, range, deceleration, 
maximum deceleration, time rate and maximum time rate of 
average and local stagnation region heat transferred into vehicle, 
total heat transferred into vehicle, and time of flight are given in 
a similar manner as those given by Allen and Eggers for great 
circle flight at small angles of inclination. 
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Symbols 


reference area for lift and drag expressions, sq. ft. 
dimensional constant in heat-transfer equations 
drag coefficient 

equivalent skin-friction coefficient Cp’ = 


1 Pl Vi Cp, 
Cr, (—)}{ — >, ( as 
s (Mr) ve 


skin-friction coefficient 
lift coefficient 
specific heat at constant pressure 
drag, lb. 
acceleration due to force of gravity, ft./sec.? 
convective heat-transfer coefficient, ft.-Ib./ft.2 sec. °R. 
convective heat transferred per unit area, ft.-lb./ft.* 
a constant = 6.8to15 X { x 
(given in references 3 and 4) 
thermal conductivity 
mass of the vehicle, slugs 
distance of glider measured from center of minor circle, 
ft. When altitude of the glider is small in compari- 
son with minor circle radius, ry may be taken approxi- 
mately as 7% 
radius shown in Figs. 1 and 2, ft. 
radius of curvature of flight path, shown in Fig. 2, ft. 
convective heat transferred, ft.-lb. 
range, ft. 
glider distance from center of earth, ft. 
radius of earth, ft. 
distance along flight path, ft. 
surface area, ft.? 
time, sec. 
temperature, °R 
velocity, ft./sec. 
Prandtl number 
air density, slugs/ft.* 
coefficient of viscosity, lb.-sec. /ft.? 
angle of inclination or angle of flight path to local 
earth horizontal, degrees, positive for descent 
surface emissivity 
angle defined in Figs. 1 and 2 
Stefan-Boltzmann constant = 3.7 X 107" ft.-lb./ft.? 
sec. °R.4 
nose or leading edge radius of body or wing, ft. 
ratio of specific heats 
angle defined in Figs. 1 and 2, degrees 
constant in density-altitude relation p = 
Here p) = reference density slugs/ft.* 6 
stant and y = altitude 
constant = (L/D)(1/2)(CpA/m)(sin \o/8) 
constant = cos p/(L/D)(1/2)(CpA/m)Ro sin? Xo 
constant = 


constant = cos 6; + a(1 — b)p,; 
constant = sin? Ao 

constant = (1 — 6)" 

constant = 2(D/L)n’ 

constant = (CpA/m)po 


pot” BY. 


= 


C = 
Cr’ = 
a 

a Cp 
D = 

h = 
St 
H = 
= 

“4 

= 

Q = 
R = 
R, = 
R = 
= 

co 
se 
T = 
P, = | for 
p = an 
6 
se’ 
av 
do = lift 
the 
ne, 
Ck 
= 
re- 
B = cir 
abl 
- 
rat 
ini 
b = 
1S 
acc 
of 
the 
Ch 
n => 
are 
GQ = 
= 
= 


cal 


C; = constant = (8/sin? Ao) [cos 6; + (a/n’)p;) 
C;' = constant = B cos 6; + aBp;(1 — b) 

Cs = constant = n’/ap 

= constant = [py + Capo cos 6, ] 
Cs = constant = (1/4)Cp’SCippV 9! 

C; = constant = (1/n’ sin \9)(1/V;) 

Cs = constant = (1/4)Cr’(p; + cos 
Cy = constant = (1/4)Cpr’CipoV;* 

Ci = constant = py + Cipp cos 6; 

Cu = constant = 

C2 = constant = Ko~'/?V,3 

F B 1 
= sin? Xo sin — cos 6] + Cig 


(preferable to use when 6 is large) 


(Go) 


8 1 | — 
= E + (C8)? + (G8) ] E= 


(preferable to use when @ is small) 


Subscripts 


f = condition at end of power boost or condition at be- 
ginning of unpowered glide 
= local conditions 


r = recovery conditions 

s = stagnation conditions 

s.l. = sea level conditions 

av = average values 

i = initial conditions 

max = maximum values 

w = wall conditions 
Introduction 


oo INTEREST on boost glide or orbital re-entry 
into planetary atmosphere has resulted in the 
consideration of three limiting factors; they are (1) 
severe vehicle deceleration, (2) severe aerodynamic 
load, and (3) intense aerodynamic heating. Besides, 
for tactical and strategic reasons, quick descent at large 
angles of inclination or forced flight at a minor circle 
path has also received considerable attention. For 
several special types of entry, analytical theories are 
available. In the case of ballistic-type entry without 
lift at sufficiently large angles of inclination where both 
the gravity force and centrifugal force are being 
neglected, the analysis of Gazley,' Allen,? Eggers,’® and 
Chapman‘ are available. In the case of gliding type 
re-entry at small angles of inclination along a great 
circle path, the analysis of Allen and Eggers’ is avail- 
able. For re-entry of a satellite with a small lift-drag 
ratio or re-entry of an orbiting vehicle with a small 
initial angle of inclination, the analysis of Chapman‘ 
is available. In the case of specific vehicles, more 
accurate machine computations® ° are available. How- 
ever, most of the analyses are limited to either the case 
of nonlifting vehicles at large angles of inclination or 
the case of lifting vehicles at small angles of inclination. 
Chapman’s solutions are much more general but they 
are primarily for great circle flights. It is the purpose 
of this paper to develop a crude approximation, but 
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Cis = constant = (4g/8V;2)n 
Ci = constant = 1 — (4g/8V;?)nF,(6;) 


1 +5 + EV(C? + 3*) — 6 zt 


Cis = constant : 


(C2 + 1) (C2 + 3?) 


C?+3 (ari) 
E = constant = (8/C; sin? \»)? = 1 
Cig = constant = (1/2) + 2?)|=> 
(1/2)[1/(Gi? + 
Cy = constant = (1/2)E[1/(C,? + 3?)] = 
+ 3?)] 


- C; sin 26 — cos 20 |<: sin? 6(C, sin @ — 3 cos ot 


240C\4 120C,4 | 12C,2 (G10) — 12064 * 1262 


for a more generalized case; and then reduce it to each 
individual simplified case for comparisons wherever 
obtainable and applicable. For great circle flight at 
small angles of inclination, the present results’ * are 
precisely reduced to those of Allen and Eggers.’ For 
nearly zero lift-drag ratio with a small initial angle of 
inclination (which becomes gradually larger as it 
enters atmosphere), the present results" '! approximate 
those of Chapman‘ and those of Rubesin-Goodwin’s 
exact machine computations. For equilibrium tem- 
perature flight, the results agree well with others.‘~® !” 
Besides these special cases, this paper presents a sum- 
mary of results on (1) minor circle flight at small angles 
of inclination,’ * (2) minor circle flight at large angles 
of inclination,'' (3) constant aerodynamic load flight, 
(4) constant equilibrium temperature flight, (5) con- 
stant vehicle deceleration flight, (6) constant vehicle 
descent flight, and (7) constant heat input flight. The 
approximate solutions developed here for various kinds 
of re-entry are approximately useful to study the de- 
celeration, heating rate and total heat absorbed for 
entry into Venus, Mars, and Jupiter, as well as for 
entry into earth atmospheres. 


Analysis of Dynamics 


Hypervelocity vehicles may either take off from 
earth by rocket boost or take off from satellite orbits 
by rocket control. These vehicles designed to reach 
their destinations will, in general, fly along a great circle 
path. However, the possibility of reaching these 
destinations by a controllable flying path other than a 
great circle path may be of interest. Recognizing 
that the great circle path is the shortest path between 
any two given points and recognizing that the great 
circle path is also the vehicle’s natural automatic glid- 
ing path (in contrast to the minor circle path which 
needs a continuous accurate aerodynamic control to 
force the vehicle flying along the desired path); never- 
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theless, one might use a minor circle path simply for 
tactical and strategic reasons. Since a minor circle path 
automatically reduces to a great circle path, when the 
minor circle angle ‘‘\y’’ is 7/2; the conventional great 
circle flight is considered here simply as a special case 
of the present more generalized case of minor circle 
flight. With this in mind, the following assumptions, 
similar to those made by Eggers and Allen* and Loh,’ 
are being made: 

(1) Flight is in planes containing the minor circle arc 
between take-off and landing (for boost glide) or be- 
tween re-entry and landing (for orbital re-entry). 

(2) Constant gravity acceleration, i.e., g= go. Note 
that this also implies the assumption that the altitude 
of the vehicle, when compared with earth radius, is 
negligible. Therefore, R, = Ro, r = ro, \ = Xo. 

(3) Rotation of earth is neglected. MINOR CIRCLE 

(4) Contribution to range of powered phase of flight FLIGHT PATH EARTH SURFACE 
is neglected. 

With these assumptions and referring to Figs. 1 and 
2, one obtains immediately the following equations of 


motion normal and parallel to the direction of flight: 


MINOR 
CIRCLE 


Fic. 1. Three-dimensional view on minor circle flight. 


L cos a = mg cos X (la) 
L sin a — mgsin \ = —m(V?/r,) (1b) 
—D + mg sin d sin 6 = m(dV/dt) = (1/2)m(dV?/ds) 
(Ie) 3 Sine 
Fic. 2. Two-dimensional view on minor circle flight. 
—1/r, = — 0)/ds], —dy/ds = cos6/r (1d) 


Combining Eqs. (la), (1b), (1c), and (1d) and using the approximations as listed, one obtains (after certain sim- 
plifications) : 


dV? ‘ 2g o(1 ) ye (“Vy 
(L/D) cos* Xo sin Ag COs he — 2g sin sin = 0 2 


The necessary aerodynamic control at any moment along the minor circle path is obtained here as a function of the 
the instantaneous velocity ‘“V,” instantaneous angle of inclination ‘‘@’’ and minor circle angle ‘‘\o”’: 


1 <6 
sin Ao cos 6 ( | 
sina gRo sin? Xo g sin Xp cos 6 \ds, 


gRo sin® Xo g sin Xo cos @ \ds 


In order to fly the the vehicle along the desired minor circle path, the automatic servo-control must be designed in 
such a way that the aerodynamics control or other control forces continuously satisfy Eq. (3). 


Analysis of Thermodynamics 
material used as heat sink. The second parameter is 


a measure or indication of the required average and 
peak average flow rate of coolants. The third param- 
eter is a measure or indication of the required local 
coolant flow rate and local structural strength. In 
developing these parameters, several crude approxi- 
mations* were used, including: 

(1) Convective heat transfer alone was considered. 
This is possible because the radiant heat transfer does 


For a generalized case study, such as the present one, 
the parameters on aerodynamic heating most suitable 
to use are those proposed by Allen and Eggers,* which 
are: 

(1) The total aerodynamic heat input to the vehicle. 

(2) The time rate and maximum time rate of average 
heat input per unit area. 

(3) The time rate and maximum time rate of local 
nose zone or stagnation zone heat input per unit area. - 

The first parameter is a measure or indication of the * For more complete discussion of these approximations, see 
required coolant weight, either liquid coolant or solid Ref. 3. 
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not appreciably influence convective heat transfer; 
and, therefore, one can be studied independent of the 
other. 

(2) Prandtl number is a constant (Cpy)/k = 
stant. 

(3) Reynolds analogy is applicable (4/pVCp) = 
(Cp/2). 

(4) The recovery temperature is much larger in com- 
parison to the wall temperature 7, — 7,=7,—-—T7= 
T, & (V2/2Cp)V P,. 

(5) Shock wave boundary-layer interaction is neg- 
lected. 

(6) Gaseous imperfection is neglected. 

Naturally, these assumptions, as explained by Eggers 
and Allen,* are not permissible for an accurate quanti- 
tative study of a specific vehicle; but they do not in- 
validate a comparative analysis which is intended to 
yield information for a comparative generalized case. 
On the basis of the foregoing assumptions, the following 
equations may be written out immediately :* 


dil Cr V? 
F = kh(T, — ( = pVCr) ( vP,) 
2 


dt 2Cp 


con- 


dQ dH | 
(4 
dt s dt halite 4) 


. Pi (=) 
C,’ = VP,(—!) ds 
(ENG Cp ‘ 


and is here assumed to be a constant.* From Eq. (4), 
the total aerodynamic heat input becomes: 


dt dt = ds = fo ds 
(5) 


By definition and Eq. (4), the time rate of average heat 
input per unit area* becomes 


dH. 1 (dQ 1 
dt 3 (2) 6) 


Here 


But more important than the average heat input per 
unit area is the maximum or stagnation point heat in- 
put per unit area. Because it seems unlikely that a 
pointed nose or leading edge will be of practical interest 
for re-entry; blunt nose alone will be considered for 
stagnation zone heat transfer. The analysis developed 
by Lees, Chapman, Romig, and at Avco*® ”~" is 
generally considered to be fairly satisfactory to use for 
re-entry of velocities up to the satellite velocity and for 
altitudes up to the extreme limit of the continuum flow 
region. Combining it with the limiting cases of heat 
transfer that occur in the free molecule flow at very 
high altitudes permits a fairly accurate calculation of 


* Eggers, Allen, and Neice* found that the effects of Mach num- 
ber and Reynolds number variation are nearly compensating. 
The variations in Cr’ for typical cases are within a reasonable 
bound. Therefore, the assumption may, in most cases of interest, 


be justified for comparative purposes. 


the heat transfer to a specific vehicle entering the at- 
mosphere. However, for the present generalized 
cases, only the theories of continuum flow, which is 
our major interest of re-entry, will be considered. For 
high-speed flight in continuum flow, the compression 
of the gas in the shock wave preceding the blunt body 
can raise its temperature to several thousand degrees 
and, thus, cause dissociation and ionization of the gas 
molecules. Under these conditions correlations for 
heat transfer at stagnation point were formulated by 
many authors. However, the accuracy of these corre- 
lations lie in the present inadequate knowledge of the 
transport properties of high-temperature gas. The 
general correlation’ ® '*—! is given as: 


dH, C ( p V ) Btu ft.-2 (7) 
= - - u ft.~? see. 7) 
dt Vo Ps.t. V 


Where the constants C, 7, and 7 depend on the type of 
boundary layer flow and transport properties of the gas. 
Since the Reynolds numbers involved during re-entry 
are usually relatively small, laminar flow is, in general, 
maintained throughout most portions of the re-entry. 
In such cases, 7 = 1/2. The value of C varies from 
16,000 to 20,000* '?~' for air and a value of 17,000 
Btu ft.*/* sec.~' for air was suggested by Chapman.* 
For gases other than air, the theory of Lees* \ gives 
C~ V — For sim- 
plicity reasons, 7 = 3 is also used here;* ‘ this value 
of j corresponds to a gas with viscosity proportional to 
7°. After putting into Eq: (7) the suggested values 
of C, i, and j, one obtains the time rate of nose zone or 
stagnation region heat input per unit area.* 4 


dH, y— 1\'4 
= K Rac ( ; ) | x 


V3 K’ V3 (8) 
o 


Following Ref. 4, it is noted that the expression in the 
bracket is a function of the planet and its atmosphere; 
for our present purposes, it may be assumed to be a 
constant for a given planet. The remaining terms 
represent the geometry of the vehicle and its particular 
type of trajectory which dictates p and V relationship 
according to the design values of (L/D), (CpA/m), and 
ao. Based on this simplification, the time rate of stag- 
nation region heat transfer per unit area may be writ- 
ten*. 4 simply as 

(dH,/dt) = K'V(p/o) V? (Sa) 


Flight at Constant Positive Lift-Drag Ratio and 
Small Angles of Inclination 


At small angles of inclination, the following approxi- 
mations! * +7 may be made: 


sin =O cos 6 = 1 (dé/ds) = 0 


Eqs. (2) and (3), therefore, reduce respectively to: 


a 
an 
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< 
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+ do + sin? A 0 (9) 
ds (L/D) gRy sin” Xo 


1 — (V?/gRo sin? Ao) 


sin a & sin 
V cos? ho + sin? Ao[l — (V?/gRo sin? Ao) |? 


(10) 


(1) Range 


Eq. (9) may be rewritten as: 


d\ (11) 
ds = — ) 
2g \D/ V1 — 2(V2/gRo) + (1/sin® do)(V2/gR,)? 


Integrating Eq. (11), one obtains immediately the expression for flight range: 


{ V2 1 V2 2 V2 
—2 1 


Ry\{L 
R=s= ( = sin Ay In (12) 
gRy Xo sin gRo sin Ao gRo 
(2) Flight Path Angle 
For exponential planetary atmosphere*:!~® p = pe” (13) 
The following equations may be written: 
—B(Rg— Ro) (B/sin Xo) (r — ro) 8 B 14 
p = pol = pot dp = - pdr = — p sin 6 ds (14) 
sin Ay sin Ao 
d 
pV sin (15) 


dt sin Ag 


Combining Eqs. (9), (10), and (15) one obtains after certain simplifications: 
Ry 
Ro @ B (: ( cot? No ; 
D gR, 


(3) Deceleration and Minimum Deceleration 


Deceleration may be obtained by rewriting Eq. (9) ina slightly different form: 


g\dt/ (L/D) gR, gRo 


Minimum deceleration may be obtained by differentiating Eq. (17) with respect to “V’’ and setting it equal to 


“0.” This results in: 
1 (dV 
= + —— (18) 
g \ dt J min cot? (1 + cot? Ao)? 


dH = =) (19) 
‘dt min 1 + cot? Xe 


(4) Altitude 
Combining Eqs. (9), (13), and (14), one obtains, after simplification : 


1 1 +] gRo gRo (20) 
eR, — gRy 


4. 


(5 


(6 


(7; 


M 
to 


(8) 


Di 


T 

" 
v2 
( 

+ 
a * This density relationship is based on the assumption of an isothermal gas in a uniform gravitational field. 63 
a 


(9) 


10) 


3) 
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|? 3 |? 2 
2 Av 
Ak 
2\D m ) 


The distance ‘“‘r”’ of the vehicle from the center of the minor circle may be obtained from the approximate relationship: 


(20’) 


l 
R, Ry => (r ry) (21) 
sin Ay 


(5) Time of Flight 


Time of flight may most easily be obtained by integrating the following expression graphically : 


gRy 


4g \D ( 3 V2\2 V2 (22) 
~ + Gr) 
sin? Ao/ \gRo gR, gR, 
Here the expression for ‘‘ds’’ is obtained from the basic Eq. (9). 
(6) Total Aerodynamic Heat Input to the Vehicle 
Combining Eqs. (11), (15) and (20), and performing integration between limits one obtains: 
m 
(23) 
( [V; 


(7) Time Rate and Maximum Time Rate of Average Heat Input per Unit Area 


Combining Eqs. (6) and (20’), one obtains: 


dH ay 1 Cp’ mgV gRy ( ( V 
~ ot? Ao 
( dt ) 2CpA (L/D) oR, T oR, co V eR, (24) 


Maximum time rate of average heat input per unit area may be obtained by differentiating Eq. (24) with respect 
to “I” and setting it equal to “0.” This results 


dH,y\ 1 C,’ mgV gR, ( (=) ( Vi ) 
a= t? Ao (29) 
( dt 2 CpA (L D) gRy V oR, 


_ \21—V1 — (3/4). F cot? Av) 
Vi = 26 
3 + cot? Ao \ ‘Vek 


(8) Time Rate and Maximum Time Rate of Nose Zone Stagnation Region Heat Input per Unit Area 


Combining Eqs. (Sa) and (20’), one obtains: 


dt gR) + gR, \ gRy 4) 


Differentiating Eq. (27) and setting to zero, one obtains: 


= K ( (1 ( ) rot? Ao ) 28 


(ae) jol (24/25)(1 + cot? Ao) VoR. (29) 
dt J mar ~ 1G 1 + cot? 
Special Case Where \y)=7/2 
When » = 2/2, the minor circle coe a great (26), (27), (28), and (29) reduce to the following simple 
circle; and Eqs. (12), (16), (17), (20), (23), (24), (25), forms, respectively : 
* When is smaller than the quantity, the maximum occurs ) In | gh) ] (12’) 
at the start of the unpowered flight. 2/\D 1 — (V//gRo) 


11) 
4) 
6) 
4 
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Fic. 3(a). Comparison of present approximate analysis with 
more exact machine calculations for ARDC model atmosphere. 
Angle of descent and altitude-velocity trajectory. 


aha: — NUMERICAL INTEGRATION FOR ARDC 
---~ CHAPMAN'S ANALYSIS 
© PRESENT ANALYSIS 
4 
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435 220 
WY 
S3 
a \ 
= 9 As \ 
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0 2 4 6 8 10 12 
DIMENSIONLESS VELOCITY 
Fic. 3(b). Comparison of present approximate analysis with 
more exact machine calculations for ARDC model atmosphere. 
Deceleration and distance traveled. 


sin = (gRo/ V2) (16’) 


(17’) 
(=) | 
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Q = (m/4)(Cr’S/CpA)(V? — V?) (23’) 
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dt 2 (L/D) 
(4) _ Cr’ _mg 25’) 
dt Jnaz 30/3 CPA (L/D) 


V = V dHav R 26 : 


dH, 2 2mg 
( dt 33 (g 0) ( ) 


2 
| if "(dHs = R 29 
( dt ) mar ( 


Eqs. (12’), (16’), (17’), (20”), (23’), (24’), (25’), (26), 
(27’), (28’), and (29’) are precisely the solutions of Allen 
and Eggers.” * 


Flight at Constant Negative Lift-Drag Ratio and 
Large Angles of Inclination 


At large angles of inclination and negative L/D, use 
is made of an assumption that the term 


cos 6{1 — (V?/gRy sin? Xo) | (30) 


is relatively small compared with other terms; and, 
consequently, it may be neglected. The justification 
of this assumption at large angles of inclination is based 
on the following reasons: 

(1) At initial phase of re-entry, the velocity ‘| is 
close to circular orbital velocity; so that 


[1 — (V2/gRy sin? X»)] & 0 


(provided Ag is relatively large, say close to 1/2, which is 
usually the case of our major interest). 

(2) At terminal phase of re-entry, the angle of in- 
clination ‘‘@’’ is close to 7/2; so that 


cos cos(r/2) = 0 


(3) During the middle portion of re-entry, ‘‘cos 
6{1 — (V?/gRpo sin? Xo) remains small, because of the 
offset effects of cos and [1 — (V?/gR, sin? Xp) |. 

Fig. 3 shows errors introduced, when L/D = 0), by 
this approximation derived only for negative L/D (in 
this case 0.1) when compared with positive L/D = 0.1. 

Based on this assumption, Eqs. (2) and (30), there- 
fore, reduce respectively to (L/D = absolute value): 


ds (L/D) ds (V?/g)(d0/ds)/ 


2g sin Ay sin@ = 0 (31) 


V2 
32 


(1) Flight Altitude and Altitude Density 


Combining Eqs. (1b), (Ic), (14), (31) and (32), 
one obtains after simplification: 


sir 
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sin 6d9 = L\(CoA\ sin ro x Combining Eqs. (13) and (35), one obtains flight alti- 

tude: 

( y = —(1/8) In + Ci(cos 6, — cos (36) 
(L/D)(CpA/2m) p( V?/g) dp (33) 


(2) Flight Velocity 
For cases of our major interest, \» is assumed to be large 
or close to 7/2, the second term in the bracket generally 
is a minor term and its magnitude usually is small in 
comparison with unity; an approximation of treating dV? 4 (5) (37) 
in this bracket as a constant according to the basic dé (L/D) 
assumption of Eq. (30) may then be made. With 
this approximation, Eq. (33) may be approximately 
integrated. The result after approximation and sim- 
plification is: 


Using Eqs. (14), (31), (35), and neglecting small 
terms, one obtains: 


Eq. (37) may be solved by the conventional method of 
solving the linear differential equation; and its solution 
after simplification becomes: 


V p? — 6? — dbtan-' (1 b)V p? — b? = (c — cos 
(34) (3) Flight Distance Along Flight Path 

Using Eqs. (lc), (31), (35), and neglecting small 


An alternative crude approximate expression of Eq. é 
PP P 4 terms as before, one obtains: 


(34) may be obtained from Eq. (33) directly by neglect- 
ing the small term in the bracket entirely: 


(39) 
sin —_ (8/sin? No) cos 6} 


p = (c’ — cos 0)/a (35) 
After integration between limits, one obtains: 
When C;” > (8/sin? Ao)? 


sin dy 39a) 


sin? 


When C;? < (@/sin? do)? 


s= In 
sin Xo (; ) £5 ( ) 6 
1 sin? Ay sin” Ay sin” A» Xo 


(40) 
(4) Flight Range 
From Figs. 1 and 2, one obtains 
dR = cos 6 ds (41) 
Sunilarly, combining Eqs. (39) and (41), one obtains after integration between limits: 
C; sin? 1 — bd) sind 
R= ( k ) sin (6 — 6,) (42) 
8 
(5) Flight Deceleration and Maximum Deceleration 
Deceleration may easily be obtained from combination of Eqs. (1b), (1c), (35) and (38): 
— E + (cos 6; — cos + + g sin ro sin (43) 
dt 2 Po sin do 


Maximum deceleration may be obtained by differentiating Eq. (43) with respect to ‘‘@’’ and setting it equal to “‘0.”’ 


This results in: 
dV V/ ( BC ) | Ci(6s—6:) 
— 0, — 6 44) 
E + (cos 6, — cos | ( 


IR 
| 


* Minor terms being neglected in the derivation. 
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C2Cy’ — — + BY 


(6) Time Rate and Maximum Time Rate of Average Heat Input per Unit Area 
Combining Eqs. (6), (35), (38), one obtains: 

LET 
) = (Cs — Cocos 8) [CisFi(8) + CeO (46) 
dH 
dt ) =* (Cy — Cy cos (47) 
6. = =* cos—! ~ 3 € s) (48) 

dt max 


‘ 


“BIC? + (2/3)2) 


(7) Time Rate and Maximum Time Rate of Local Stagnation Region Heat Input per Unit Area 


Combining Eqs. (Sa), (35), (88), one obtains: 
dH, 


it = Ci(Cw — Cu cos (8) + Cag (49) 
at 
(=) >* Ci2(Cio Cu cos 6;) e(3/2)OF— (50) 
max 
— 01") + (5.8) 
6; = cos! 3 3 (51) 
dt } max — 


BIC? + (1/3)?] 


(8) Total Aerodynamic Heat Input to the Vehicle from the Beginning to the Point Concerned 


Combining Eqs. (5), (14), (35), (38), one obtains: 


Gs Cs cos 0 


(° 
Q =*(—— dg 
sin Ag of B (52) 
sin® 


(9) Flight Time 
| 16 6 (Ci/2) 16 
t= fe = aj [CisF\(0) Cie — 
os 


Here the expressions of ‘“‘ds’’ and “VV” are obtained 
from Eqs. (38) and (39) respectively and the last two 
expressions can easily be integrated by graphical means. 
For high-speed re-entry, except the terminal phase of 
flight, the gravity term may sometimes be neglected. 
In such cases, Cj; and C\4 become 0 and 1.respectively; 
and the foregoing equations are all reduced to much 
simpler forms. 


Special Case Where )\y = 7/2 


When XA» = (2/2), the minor circle becomes a great 
circle; and all the results obtained here are simplified 


* Minor terms being neglected in the derivation. 


C; ( 2 -) cos @ C; ( ) cos 6 
sin” Ao sin? Xo 


and are reduced to the corresponding results valid for 
great circle flights.'’ Since no analytical or numerical 
results were found in the literature even for comparing 
order of magnitude for large angles of inclination at 
constant negative lift-drag ratio, a case of nearly zero 
lift-drag ratio was used here for this purpose where the 
term lift-drag ratio, whether positive or negative, is so 
close to zero that the trajectory itself is very close to the 
trajectory of zero lift-drag ratio. Fig. 3 shows such a 
comparison and gives errors introduced, when L/D = 0, 
from the present crude approximation at negative lift- 
drag ratio of 0.1 and those more accurate and exact 
solutions for positive lift-drag ratio of 0.1 by Chapman 
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and Rubesin-Goodwin. The differences were small as 
shown for the two cases, as might be expected for 
L/D = 0 (except deceleration). 


Trajectories of Variable Lift-Drag Ratios* 


Constant lift-drag ratio trajectories were discussed in 
the preceding sections. Constant lift-drag ratio re- 
quires variable angles of inclination. In this section, 
trajectories of variable lift-drag ratios for lower atmos- 
phere penetration will be presented. Variable lift- 
drag ratio may be obtained, after a certain phase of 
penetration into lower atmosphere, either by variable 
angles of attack or by use of a body the drag coefficient 
of which remains essentially constant when its lift 
coefficient varies (bodies with high parasite drag co- 
efficients which overdominate induced drag _ coeffi- 
cients). However, no matter what kind of scheme is 
used, the servomechanism must be designed in such a 
way so that the required variation of lift-drag ratios 
(by aerodynamic or reaction) is assumed to be always 
satisfied during the entire specified lower atmosphere 
trajectory. Trajectories of major interest include the 
following: 


(1) Constant deceleration flight 

(2) Constant aerodynamic load factor flight 
(3) Constant rate of descent fight 

(4) Constant angles of inclination flight 
(5) Constant rate of heat input flight 

(6) Constant equilibrium temperature flight 


For simplicity reasons, only the great circle path will be 
considered. For great circle path, \y = (2/2), cos Ay = 
0, sin Ay = 1, Eqs. (2) or (1b) and (1c) are considerably 


simplified and they are reduced to:7 
dV? 4 | (: V? ) Vy? 9 
C == 
ds (L/D) g \ds 


or 


_dV_ idl (55) 
m dt 2ds 
L (1 ) 5 
= —mvV* 56 
mg COs (56) 


Flight at Variable Angles of Inclination 


Based on the foregoing Eqs. (54) or (55)—(56), to- 
gether with equations (5), (6), (8), the solutions on in- 
stantaneous velocity, density, distance along flight 


* Such flights are usually considered for a portion of re-entry 
trajectory for the purpose of reducing aerodynamic heating or 
deceleration or other purposes. 

+ For high-speed re-entry, the gravity component in the drag 
direction is usually very small in comparison with drag itself; 
consequently, it may be neglected.'~* This does not apply at 
the terminal phase, or the case where the gravity component is 
important. 


path, range, deceleration, maximum deceleration, time 
rate and maximum time rate of average and local stag- 
nation heat transferred into vehicle, total heat trans- 
ferred into vehicle and time of flight may be solved 
similar in manner to those presented previously. For 
simplicity, only the final results are summarized in 
Table 1. Derivations are self-explanatory there. 


Flight at Constant Angles of Inclination 


At constant angles of inclination, Eq. (54) further 
simplifies to: 


= 57 
ds (L, D) 


or 


V2 
L = mg cos 6 (1 _ =) (58) 


dV idv?.. 
V? = => 
( m dt 2 ds 


The solutions on instantaneous velocity, density, alti- 
tude, distance along flight path, range, deceleration, 
maximum deceleration, time rate and maximum time 
rate of average and local stagnation region heat trans- 
ferred into vehicle, total heat transferred into vehicle, 
and time of flight are also presented in Table 1. A 
typical case of the present solution is for (L/D) = 0. 
When (L/D) = 0, the vehicle becomes a nonlifting 
body, such as a nonlifting missile. The present solu- 
tions then yield the solutions obtained by Gazley,' 
Allen and Eggers,” and Chapman‘ for nonlifting bodies 
entering earth atmosphere at large angles of inclination. 


Discussions 


Because of the space limitation of this paper, only 
significant points of re-entry dynamics and thermo- 
dynamics will be mentioned briefly.'~' 


(I) Small Angles of Inclination vs. Large Angles of Inclina- 
tion at Constant Lift Drag Ratio 


(1) In general, small angles of inclination will result: 


(a) Deceleration accomplishment at higher altitudes 
thus entering lower atmosphere at lower velocities. 

(b) Lower heating rates. Heat transferred into 
vehicle per unit time is less. 

(c) Smaller deceleration. Time for deceleration is 


longer. 

(d) More heat absorbed. Total heat absorbed into 
vehicle is increased, because a larger value of C,’ is 
resulted. 

(e) Generally positive lift drag ratio required for re- 
entry at small angles of inclination. 


t CorreEcTION: In Table 1, the algebraic sign of the first 
term on the right-hand side of the (L/D) expression, as shown in 
the second, third, fourth, and fifth columns, should be reversed. 
The second term on the right-hand side of the (L/D) expression 
as shown in the fifth column of Table 1 should be multiplied by 
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Trajectories of Variable Lift-Drag Ratios 


(Used primarily for lower atmosphere penetration for purpose of modulation of aerodynamic heating, deceleration, etc.) 


FLIGHT AT CONSTANT DECELERATION 


FLIGHT AT CONSTANT RATE OF AVERAGE 


TRAJECTORIES OR FLIGHT AT CONSTANT AERODYNAMIC HEAT INPUT PER UNIT AREA OR CONSTANT 
LOAD FACTOR AVERAGE RADIATION EQUILIBRIUM SKIN 
Constant =K,=- dH 3 4 
= Constant f; Wave Vv; 
-(- 2m Vi yy 
vetocity \/ \/=( Fora) e Dy Ve eff 


ALTITUDE Y 


ALTITUDE DENSITY 


FLIGHT DISTANCE ALONG 
FLIGHT PATH yg 


RANGE “R” 


via 


dv 
DECELERATIO! ( 


dt 


MAXIMUM DECELERATION 


(FE) mas: 


MAXIMUM DOES NOT EXIST 


MAXIMUM DOES NOT EXIST 


VELOCITY AT 
MAXIMUM DECELERATION 


DOES NOT EXIST 


DOES NOT EXIST 


RATE OF AVERAGE HEAT ) 
IMPUT PER UNIT AREA dt 


KY 


=k, 
dt 


MAX. RATE OF 


dH 
AVERAGE HEAT ( at max, 


IMPUT PER UNIT AREA 


MAXIMUM DOES NOT EXIST 


MAXIMUM DOES NOT EXIST 


VELOCITY AT MAX. RATE OF 
HEAT IMPUT PER UNIT AREA 


DOES NOT EXIST 


DOES NOT EXIST 


RATE OF STAGNATION 
POINT HEAT IMPUT (2) 
PER UNIT AREA 


MAX. RATE OF STAGNATION 
POINT HEAT IMPUT ELA 


PER UNIT AREA t /max. 


MAXIMUM DOES NOT EXIST 


MAXIMUM DOES NOT EXIST 


VELOCITY AT MAX. RATE OF 
STAGNATION POINT HEAT IMPUT 
PER UNIT AREA 


DOES NOT EXIST 


DOES NOT EXIST 


TOTAL AERODYNAMIC HEAT 
ABSORBED “a 


TIME OF FLIGHT 


toa 


t: 


ANGLE OF INCLINATION. Sino (Sr Sino 3 
=) $)= 3Ja [ 


LIFT DRAG RATIO (+) 


ae 


CERTAIN CONSTANTS USED 


b=4/ mC, 


THE TRAJECTORY OF CONSTANT DECELERATION 
FLIGHT IS EXACTLY SAME AS THE TRAJECTORY OF 
CONSTANT DYNAMIC PRESSURE FLIGHT (OR USUALLY 
CALLED CONSTANT AERODYNAMIC LOAD FACTOR 
FLIGHT) 


THE TRAJECTORY OF A CONSTANT AVERAGE HEAT 
IMPUT FLIGHT IS EXACTLY THE SAME AS THE 
TRAJECTORY OF A CONSTANT HEAT TRANSFER 
FLIGHT WHICH MAINTAINS A CONSTANT AVERAGE 
WALL TEMPERATURE BY EQUATING RADIATIVE 
COOLING RATE TO HEAT IMPUT RATE. 
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TABLE 1 (ConrT.) 
Trajectories of Variable Lift-Drag Ratios 


(Used primarily for lower atmosphere penetration for purpose of modulation of aerodynamic heating, deceleration, etc. ) 


FLIGHT AT CONSTANT RATE OF STAGNATION 
POINT HEAT IMPUT PER UNIT AREA OR 
CONSTANT STAGNATION RADIATION EQUILIBRIUM 
TEMPERATURE FLIGHT 


V°%= Constant v;° 


FLIGHT AT CONSTANT RATE OF DESCENT 


(2) =Constant= Ks= (22), 


FLIGHT AT CONSTANT ANGLES OF INCLINATION 


=Constant= 


V-V. ) 


cov) 


MAXIMUM DOES NOT EXIST 


avy 
dt Jmax.4 \ 7m 


DOES NOT EXIST 


Vogt max’? (Ee) 


- (2) 


(Se) Sates) 


(Fe) 


= 


MAXIMUM DOES NOT EXIS, 


‘max. 


DOES NOT EXIST 


V nas ; 


ele#)-G) 


max. 


(2 


MAXIMUM DOES NOT EXIST 


DOES NOT EXIST 


max. 


Sin 


CigCoA\ 


Sin® =- 


Sin © = Constant=Six 0; 


(FD 


a=9 mBK? 
b=( 


* 
Co ) 


V; 


THE TRAJECTORY OF CONSTANT STAGNATION POINT 
HEAT IMPUT FLIGHT IS EXACTLY THE SAME AS 
THE TRAJECTORY OF A CONSTANT HEAT TRANSFER 
FLIGHT WHICH MAINTAINS A CONSTANT SKIN 
TEMPERATURE AT THE STAGNATION OR CRITICAL 
POINT OF THE VEHICLE BY EQUATING RADIATIVE 
COOLING RATE TO HEAT IMPUT RATE. 


A SPECIAL CASE OF THE PRESENT SOLUTION IS FOR 
WHEN S=°, THE VEHICLE BECOMES A NON- 
LIFTING BODY, SUCH AS A NON-LIFTING MISSILE. 
THE PRESENT SOLUTIONS THEN YIELD PRECISELY 
THE SOLUTIONS OBTAINED BY GRAZLEY, ALLEN, 
AGGERS, CHAPMAN (1, 2, 4) FOR WNON-LIFTING 
BODIES ENTERING EARTH ATMOSPHERE AT LARGE 


ANGLES OF INCLINATION. 
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(2) In general, large angles of inclination will result: 


(a) Deceleration accomplished at lower dense at- 
mosphere but passing it at a shorter time. 

(b) Higher heating rates. Heat transferred into 
vehicle per unit time is more. 

(c) Larger deceleration. Time for deceleration is 
shorter. 

(d) Less heat absorbed. Total heat absorbed into 
vehicle is decreased, because a smaller value of C,’ is 
resulted. 

(e) Relatively smaller lift-drag ratio or generally 
negative lift-drag ratio required for re-entry at large 
angles of inclination. 


Therefore, for minimizing total heat input or total heat 
absorbed by the vehicle for the heat-sink type re-entry 
vehicles (where radiation heat loss is small), large angle 
of inclination is required generally. But for vehicles 
operating under radiation equilibrium or for vehicles 
which have significant radiation losses, small angle of 
inclination is, in general, advantageous. Such vehicle 
can decelerate at higher altitude at smaller angles of 
inclination so that the heat input rate can be balanced 
or compensated by the heat radiated from the vehicle. 
Large angles of inclination requires to use a retro- 
rocket. When decelerations are tolerable, large angles 
of inclination are practical. Entry into Mars can 
easily be accomplished with large angles of inclination. 
Entry into earth and Venus may also be accomplished 
by large angles of inclination. With present human 
tolerances entry into Jupiter is almost impossible to 
use large angles of inclination. These are shown in the 
respective figures attached* (Figs. 5 through 12 for 
negative L/D = 0.1). 

It is to be noticed that a number of assumptions 
which were originally used for small angles of inclina- 
tions were also used here for large angles of inclination. 
Certainly, these assumptions in certain cases may not 
be completely justified in the case of large angles of 
inclination; however, for simplicity and consistency 
they were simply assumed to be also valid for large 
angles of inclination in the present generalized analysis. 
It is also to be pointed out that the present solution for 
large angles of inclination has a singular point at 
zero lift-drag ratio. Therefore it is to be expected 
that the solution probably will yield a better result 
at an absolute value of lift-drag ratio greater than 0.1. 
Finally it must be emphasized here again that the 
L/D values shown in the analysis for both small 
and large angles of inclination are their respective 
absolute values; therefore, no negative sign is re- 
quired for the L/D value when it is used for descent 
at large angles of inclination. 


(II) Trajectories of Variable Lift-Drag Ratios 


Constant lift-drag ratio trajectories were discussed 
in the previous section. Constant lift-drag ratio re- 


* By neglecting small terms, the initial phase of some skips at 
large angles of inclination are shown by its smoothed averaged 
trajectories in the present approximate analysis. 
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quires variable angles of inclination (either large or 
small). As discussed in the previous section, large 
angles of inclination result in higher rate of heat trans- 
fer but lower amount of total heat input to the vehicle; 
while small angles of inclination result in lower rate of 
heat transfer but higher amount of total heat input to 
the vehicle. Both rate of heat transfer and amount of 
total heat input to the vehicle might raise the vehicle 
temperature above tolerable limit during certain phases 
of planetary atmosphere re-entry. In such cases, 
variable lift-drag ratio trajectories may be used. Since 
The critical portion of atmosphere re-entry usually 
occurs either at relatively lower atmosphere or rela- 
tively shorter length of time, variable lift-drag ratio 
flight, although difficult, is not impossible. Special 
trajectories, after entering lower atmosphere, may be 
accomplished by variation of lift coefficient or lift-drag 
ratio aerodynamically or reactionally for a short while. 
Control of vehicle trajectories in such a way may be 
made for (1) constant deceleration flight or constant 
aerodynamic load factor flight, (2) constant rate of 
average heat input flight or constant average radiation 
equilibrium temperature flight, (3) constant rate of 
stagnation point heat input flight or constant stagna- 
tion region radiation equilibrium temperature flight, 
(4) constant rate of descent flight, (5) constant angle of 
inclination flight. Trajectories under variable lift- 
drag ratio may be controlled in such a way that the 
peak deceleration, heating, etc., experienced by the 
vehicle is held to a tolerable amount within a short 
length of time or within a certain phase of lower at- 
mosphere penetration. The summary results of the 
variable lift-drag ratio trajectories are given in Table 1. 
The following few significant results may be mentioned 
briefly. 

(1) Constant Angle of Inclination Flight—For con- 
stant angle of inclination at small lift-drag ratios, 
Ferri’s exact numerical analysis’ shows that a simple 
relationship exists 


(1, £2) (dV, => —170 sin 6 


which is independent of (CpA/m). The present closed 
form analytical result given in Table 1 reduces to 


(1/g)(dV/dt) max —174 sin 0 


In view of the numerical result obtained by Ferri and 
analytical result derived here, the agreement is con- 
sidered to be good. 

(2) Constant Rate of Heat Input or Constant Radiation 
Equilibrium Flight—For a given surface temperature, 
rate of heat radiated out is equal to 


ek’T.,' 


Using this quantity and the rate of heat input expres- 
sion, a velocity-density or velocity-altitude relation- 
ship can be determined such that the heat input to the 
body is exactly equal to the heat radiated out. How- 
ever, it is necessary to exert lift, aerodynamic or reac- 
tion, to traverse such trajectory as required for all 
the variable lift-drag ratio flights. Two radiation 
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LEGEND 
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Fic. 4. Maximum radiation-equilibrium temperature at lami- 
nar stagnation point for entry from decaying orbits into earth 
atmosphere. 
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equilibrium flights are given here; one is based on rate 
of average region heat input and the other one is based 
on rate of stagnation region heat input. Should Lees’ 
results'* be more desirable to use for hypersonic heat 
transfer, the average region rate of heat transfer may be 
obtained as a certain percentage of stagnation region 
rate of heat transfer; then constant stagnation region 
radiation equilibrium temperature flight as presented 
here may also be interpreted as constant average region 
radiation equilibrium temperature flight. Numerical 
values on constant radiation equilibrium temperature 
flight agrees fairly well with the other numerical 
results.‘~* '? This is shown in Fig. 4. 

(3) Constant Deceleration Flight—Constant decelera- 
tion results in a constant dynamic pressure ‘‘(1/2)pV,”’ 
consequently, a constant aerodynamic load factor for 
this kind of flight. This trajectory may be especially 
considered for proper entry into Jupiter where entry 
deceleration is most critical. 
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Creep Deformation and Stresses in 
Pressurized, Long Cylindrical Shells’ 


M. P. BIENIEK* ano A. M. FREUDENTHAL** 


Columbia University 


Summary 


The state of deformation and stress in a long, circular, cylindri- 
cal shell with freely supported or clamped ends, loaded by uni- 
form lateral pressure, is analyzed under condition of slow sta- 
tionary (second-stage) creep. In this analysis, use is made of an 
approximate method based on a general extremum principle for 
dissipative media. The theoretical results are checked by ex- 


periment. 
Symbols 
a = shell radius 
h = wall thickness 
l = length 
n = material constant (for creep function) 
p = lateral pressure 
w = radial displacement of the middle surface 
Wo = membrane radial displacement of the middle 
surface 

Cizy Sif = deviatoric components of strain and stress 
t; = surface tractions 
vi = velocity vector 
B,, = extensional rigidity of the shell in creep analysis 
= material constant 
D, = bending rigidity of the shell in creep analysis 
E = Young’s modulus 
j,’ = second invariant of strain-rate deviation 
J! = second invariant of the stress deviation 

q = material constant 
M, = bending moment in x plane 
Ne = circumferentia! normal force 
W = rate of energy dissipation 
X; = body forces 
a = parameter characterizing the deflection of the 

shell in creep 

B = similar parameter for elastic deflection 
Ez, & = axial and circumferential strain components 
K = curvature of the deformed shell 
A = a parameter 
v = Poisson’s ratio 
¢, v, 9, € = functions defined by Eggs. (3.2), (3.4), (3.5), (3.6) 
Or, Ty = axial and circumferential stress components 
é = dimensionless coordinate 
x = dissipation potential 
= abbreviation for a/8 
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(1) Introduction and Basic Equations 


8 bes BENDING MOMENTS arising in thin shells as a 
result of boundary restraints and of loads that 
cannot be carried by membrane stresses are of consid- 
erable importance in the design of such shells. Elastic 
theory has shown that these moments which may 
reach high intensities in the vicinity of the disturb- 
ance, display the characteristic distribution of a 
damped wave with relatively rapid decay from the 
peak intensity. Thus, it is generally assumed that only 
a very narrow strip of the shell close to the disturbance 
is affected by such moments. 

The analogy between linear elastic and linear viscous 
behavior implies that a similar distribution of moments 
will arise in a linear viscous “‘creeping’’ shell. How- 
ever, the known nonlinearity (with respect to stress) 
of creep in high-temperature metals suggests that a 
considerable reduction of the intensity of the bending 
moments, as well as a significant change in their dis- 
tribution will take place, and that this effect will be 
the more pronounced the higher the nonlinearity. 

It is the purpose of the present investigation to 
study this effect of the nonlinearity of creep on the 
bending moments in thin shells by analyzing the 
specific problem of the end moments in a thin cylin- 
drical shell under lateral pressure, since this is one of 
the problems which most frequently arises in various 
phases of high-temperature design. 

Considering a long, thin-walled circular cylindrical 
shell of uniform wall thickness /, radius of the middle 
surface a, and length /, subjected to uniform lateral 
pressure p, the only velocity component is in the radial 
direction w = w(x); the longitudinal coordinate x is 
measured from the edge of the cylinder (Fig. 1), which 
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may be freely supported or clamped. The rate of 


curvature of the axially deformed shell is 
k= Ox? ( 1) 


Introducing the Navier assumption, the longitudinal 
strain rate is of the simple form 


= 2k = 2(0°w/Ox") (1.2) 


where z denotes the (radial) distance from the middle 
plane; the uniform circumferential strain rate is 


= w/a (1.3) 


The longitudinal bending moment J/, and circumfer- 
ential normal force NV, expressed in terms of bending 
(longitudinal) stresses o, and circumferential stresses 


o, are 
h/2 h/2 
M, = o,2dz, N, = f o,dz (1.4) 


—h/2 —h/2 


Moreover, a circumferential bending moment Mg 
exists; it is related to /,, but it is not significant in the 
considered problem. 

The material is assumed to be incompressible and 
nonlinear viscous, with negligible stage of primary 
creep. In the case of uniaxial stress, the relation be- 
tween strain rate and stress is introduced, 


é = Ko" (1.5) 


where A and mw are material constants. 

In the three-dimensional case, the general relation 
between deviatoric components of the strain rate and 
stress is! 


= 


which, with f(J2’) = CJ.’~-"/? in accordance with 
Eq. (1.5), and with o, = 0, €,, = 0, may be written in 
the form 


é, = K(o,? + o,? — o20,)%—"!? (a, — 0.50,) (1.6) 


é, = A(e,* + ¢,? — (¢, — 0.5¢,) (1.7) 


Since the use of relations (1.6) and (1.7) would cause 
considerable computational difficulties, the influence 
of bending stresses on the circumferential deforma- 
tion is neglected, and the following approximate rela- 
tions are therefore introduced for circumferential 
stress and strain by setting o, = 0 in Eq. (1.7): 


é, = Ke," of o, = K-"%,"" (1.8) 


which are identical with the relations of uniaxial state 
of stress. Neglecting also the influence of the cir- 
cumferential extension on the bending deformation; 
i.e., deriving the corresponding relation under the as- 
sumption €, = 0 which implies, on the basis of Eq. 
(1.7), o, = 0.5 o,; the following relations are obtained 
for bending stress and longitudinal (bending) strain: 


(1.9) 


= = K,o," or 


For n = 1, the relations (1.8) and (1.9) are identical 
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with the relations of linear viscous flow or analogous to 
the relations of linear elastic deformation with vy = 0.5. 

If Eqs. (1.8) and (1.9) are introduced into Eq. (1.4), 
the following relations between « and /, and between 
é, and N, are obtained: 


h/2 
M, = x, edz = 


—h/2 


h/2 
f dg = Dyk!" (1.10) 
—h/2 
h/2 
where git" dg 
—h/2 


h/2 
and N, = Kom f dz = B,é,"" (1.11) 


—h/2 


h/2 
B, = K-1 f dz 
—h/2 


In the cases of solid walls 


2 
2n + 1 2 


where 


(1.12) 


(2) Extremum Principle for Stationary Creep 


The actual deformation rates in slow stationary flow 
minimize the expression? 


- tvdA XwadV_ (2.1) 
V A V 


where x(e;;) is the dissipation potential (per unit vol- 
ume), expressed as a function of the state of flow, t,v; is 
the rate of work of the surface tractions ¢;, and X jv; is 
the rate of work of the body forces X,. In the case of 
a stress—strain-rate relation in the form of a power 
function 


= Dies, or = 22) 


the specific dissipation potential x(e;;) is related to the 

specific rate of energy dissipation s;,¢;; [see reference 3]: 

n 2n 

x (éi) n+1 (2.3) 
OW _ dx(eu) _ 


so that 
n + 1 06 


as specified in Eg. (2.2). 

For the problem of the cylindrical shell, the dissipa- 
tion potential x per unit area of the middle surface 
will be used. It can be expressed in terms of /,, 
N, and &, é, in the form 


n+ 1 


+1 
or in terms of deformation rates by using Eqs. (1.10) 
and (1.11) 


(M,i+N,€,) (2-4) 


nN 


n+1 


x= (D, + Bn (2.5) 
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the rate of work of the surface forces 


t,v; = pw (2.6) 


while the body forces X¥; = 0. Hence, the expression 
(2.1) to be minimized takes the form 


n 

E " + — prs Jax min. 

2.7) 


Substituting Eqs. (1.1) and (1.3) for curvature and 
circumferential strain, respectively, the extremum 
principle which has to be satisfied by the velocity field 
is obtained in the following final form: 


9 fn D 
J, in +1 | "\ Ox? + 


wW 1+1/n” 
B, (*) | pie dx min. 


(2.8) 


(3) Analysis of the Shell Deflection 


In order to minimize the integral (2.8), some velocity 
field w(x) of the deformed walls of the shell must be 
assumed, with some parameters that can be determined 
with the aid of the minimum condition. Some pre- 
liminary experiments suggested that the deformed shape 
of the shell in creep might be fairly well represented by 
expressions similar to those obtained for the elastically 
deformed shell of the same geometry, provided the 
maximum ‘‘membrane’”’ deflection and the parameter 
a governing the wave length of the deformed generatrix 
are properly determined. This parameter for an elastic 
shell is of the form‘ 


B= V3(1 — v*)/(ah)? (3.1) 


Thus, the assumed shape of a long shell with clamped 
edge at x = 0, obtained by replacing 8 by a, is 
w(x) = wo[l1 — e-@*(cos ax + sin ax)| = 

— ¢) (3.2) 


where w, denotes the (uniform) rate of radial expansion 
at large distances from the edge. Hence, according 


to Eq. (1.11), with V, = pa, 
Wy = = a(pa, B,)" (3.3) 


The rate of curvature 


da 


With the new variable £ = ax, Eq. (3.7) becomes 


Pa) al n 1+1/n w\itin dé 
da Jo +8, (2) |- 


net nes nes 

Fic. 2. Moment and normal force diagrams for elastic shell 
(n = 1) and two cases of creep (nm = 3 and m = 5). 

x 

in 

3.0 
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1.0 

0.2 w in 

Fic. 3. Radial displacements, theoretical form = 1,” = 38, 


n = 5, and measured (broken line). 


kK = = ax — sin ax) = 
(3.4) 
Similarly, for a simply supported edge 
w(x) = — cos ax) = — 6) (3.5) 
The rate of curvature 


(3.6) 


where wy is the same as before. 


sin ax = 


i= 


the functions g(x), ¥(x), 0(x), and ¢(x) contain the un- 
known parameter a, which has to be determined from 
the condition (2.8). Consequently, 


1+1/" 
B,, (*) | dx = 0 


(3.7) 


(3.8) 


It is expedient to split the range of integration into two parts, one 0 < & < A, over which the bending moments are 
significant, the other \ < & < al, over which the work of the bending moments can be neglected. Then, with 
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al sy\ 1+1/" n 
a a a a 


Eq. (3.8) may be written in the form 


n Wo 1+1/n 
(/ d/a) B, (=) = 0 
or, after differentiation with respect to a, 
D,, 2/n d. 5 —2 peal —2 


5 1+1/n 
E Bs pi =0 (3.9) 


Making use of the assumed forms of the function w according to Eqs. (3.2) and (3.5), and considering that the de- 
flection w is determined by Eq. (3.3), relatively simple expressions for a are obtained: (a) in the case of the clamped 
edge, 


D,( + 2/n)2a dé 
0 
and (b) in the case of the simply supported edge, 
B ql + in) 6dé +f — — 1 |dé 
D,(1 + 2/n)2a ge 
0 


If the walls of the cylinder are solid and of a homogeneous material, B, and D,, may easily be determined and the 
equations for a become (a) 


(3.12) 
+ 2/01 + amy) a 
and (b) 
a= (3.13) 
+ 2/1 + ag 
which may be written in the general form 
= B&(n) (3.14) 
where the coefficient &(7) depends on the exponent 
and the type of boundary conditions at x = 0; for » = 1, w = Kp"a'+"h-"(1 — 6) 
a = B,and@*= 1. 
The final formulas for deformations, circumferential N, = pa(i - 6)" (3.16) 
forces, and bending moments for the cylinder with p 3atman — 
clamped edge can now be written M, = 38? 2 + 1/m &(n)?/"erln 
w = Kp"a't"h-"(1 — ¢) 
N, = pa(l — ¢)'" (4) Conclusions 
3.15 
The coefficient has been computed for a cylinder 
M, = 2822+ 1 /n with changed end, considering two values of For 


n = 3, © = 0.6598, and for 2 = 5, 6 = 0.5086. For 
For the cylinder with simply supported edge (Continued on page 778) 
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On the Prediction of Mixed 


Subsonic/Supersonic Pressure Distributions 


COLIN S. SINNOTT* 
National Physical Laboratory, Teddington, England 


Summary 


High-speed wind-tunnel results are analyzed to derive a semi- 
empirical scheme for the prediction of transonic pressure dis- 
tributions. The supersonic and subsonic parts of the flow are 
treated separately, and then linked by an empirical shock pres- 
sure rise relation. The significance of the empirical results is 
considered in relation to the physical mechanism of transonic 
flows. It is also shown that theoretical solutions can be im- 
proved by introducing the empirical shock relation. 


Symbols 

( = airfoil chord 

H = free-stream stagnation pressure 

M, \J, = local and free-stream Mach number, respectively 

M* = crest critical Mach number 

p = local static pressure 

x = ordinate in free-stream direction, origin at airfoil 
leading edge 

= (x — xer)/(Xs — XcR) 

a = airfoil incidence 

Subscripts 

1,2 = values immediately upstream and downstream of a 
shock wave, respectively 

SR = equivalent sonic-range value 

CR = crest value 

a = shock-wave value 

G = value calculated by Glauert compressibility rule 


(1) Introduction 


T HEORETICAL METHODS for the estimation of the 
two-dimensional inviscid flow about airfoil sec- 
tions have existed for many years. The incompres- 
sible flow about simple profiles can be calculated ex- 
actly, and useful approximate solutions for the com- 
pressible subsonic flow about geometrically defined 
profiles have been developed. However, in these 
compressible flow methods, a solution is obtained only 
by recourse to approximations to the governing equa- 
tion which destroy the feature that is essential for the 
representation of mixed subsonic/supersonic flows— 
i.e., the change from an elliptic form to a hyperbolic 
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form of the equation when the local speed passes from 
below to above the local speed of sound. This change 
in form is associated with a nonlinear term in the equa- 
tion, which, if retained, leads to great difficulties in the 
derivation of an analytical solution to even the most 
simple transonic flows; in this paper an alternative 
approach is suggested. This approach is based on ap- 
proximations that appear to be sensible from the 
physics of the flow, rather than to facilitate a mathe- 
matical analysis of the flow equations. 

Two distinct approaches to the solution of the non- 
linear partial differential equation for the transonic 
flow about airfoils have been adopted. Attempts have 
been made to calculate solutions of the exact equation 
by using the methods of numerical analysis, notably 
relaxation. References 1 and 2 are examples of this 
approach, neither of which appear wholly satisfactory. 
Emmons’ earlier work* (1946) has, however, been ex- 
tremely valuable in throwing light on the flow in the 
neighborhood of a normal shock wave adjacent to a 
solid boundary. A serious drawback in attempts to 
derive solutions to the exact equation is the amount of 
computation required. It is possible that useful re- 
sults could be obtained if electronic computers can be 
used to perform the relaxation process. The alterna- 
tive approach to the calculation of transonic flows is 
that based on approximations to the equation to facili- 
tate the use of analytical methods. It is essential in 
this approach that the mixed elliptic/hyperbolic nature 
of the equation be retained. Approximations intro- 
duced to linearize the equation, or to reduce it to a set 
of linear equations do not give a proper representation 
of the significant differences in stream-tube behavior 
in subsonic and supersonic flow. The well-known 
methods of expressing the velocity potential as an ex- 
pansion in Mach number or airfoil thickness are of this 
kind, and although they lead to solutions which include 
local supersonic regions, they are physically unreal. 
The same criticism applies to the method of solution of 
the transonic flow equations suggested by Meksyn‘ 
(1953) (see reference 5); it is interesting to note that 
this method is shown to be equivalent to the M? ex- 
pansion method by Imai® (1957). 

An analytical solution of a simplified form of the 
transonic flow equation in which the essential quadratic 
nature of the equation is retained was developed by 
Spreiter and Alksne’ (1955); the analysis follows that 
made by Oswatitsch® (1950) and Gullstrand® (1952). 
The differential equation governing the flow is con- 
verted into an integral equation to which an iterative 
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Fic. 1. Typical family of transonic pressure distributions. 


solution is obtained. Solutions for mixed subsonic/ 
supersonic flows including shock waves can be ob- 
tained, but approximations introduced to simplify 
the analysis and computation still prevent an exact 
representation of the physics, particularly in the neigh- 
borhood of a shock wave. However, Spreiter’s method 
appears to be by far the most promising theoretical 
approach to transonie flow problems, and the analysis 
of results obtained by it is described in Section (4). 
The flow considered here is that for which the free- 
stream Mach number is such that the flow over part 
of the airfoil surface is supersonic. When the maxi- 
mum local Mach number exceeds about 1.1, it appears 
that the downstream boundary conditions can only be 
satisfied in a real viscous flow by the termination of the 
supersonic region by an abrupt pressure rise to sub- 
sonic conditions. Although this type of flow is cer- 
tainly a possible theoretical inviscid flow, no entirely 
satisfactory explanation of its predominance over the 
sometimes-possible smooth supersonic deceleration has 
yet been found; a number of hypotheses are discussed 
in reference 10. Nevertheless, in this paper, it is sup- 
posed that all limited supersonic regions terminate in 
shock waves. The shock-wave position must then be 
such that the pressure on the airfoil surface immedi- 
ately upstream of the shock, (/,), and that immedi- 
ately downstream, (/2), which must be related via the 
flow over the rear of the airfoil and the wake to the 
downstream boundary conditions, are such that the 
ratio p/p, satisfies a consistent relation. Analysis of 
the stream tube adjacent to the surface suggests that 
the shock pressure rise on the surface should satisfy the 
classical normal shock pressure rise relation. Now, 
in addition to the complexity of the transonic flow equa- 
tion mentioned earlier, a further problem in the corre- 
lation of theoretical and experimental studies of trans- 
onic flow is the finding that measured shock pressure 
rises are invariably less than those given by the shock 
equation. Moreover, they are not found to follow 


any definite trend. Fortunately, the detailed numeri- 
cal solutions for mixed flows in a hyperbolic nozzle ob- 
tained by Emmons’ suggest an explanation of the dis- 
crepancies. Emmons found that if he used the normal 
shock equation to determine the pressure change 
through a shock, the predicted pressure rise was fol- 
lowed immediately by a pressure fall. Measurements 
of the pressure changes through shock waves made by 
Ackeret, Feldman, and Rott!! show a similar vari- 
ation. Emmons explains this behavior by considering 
the effect on the pressure gradient normal to the sur- 
face of applying the shock relation along a line normal 
to the streamlines. It is clear that a reversal in direc- 
tion of the normal pressure gradient must occur, and, 
hence, there must be a discontinuity in surface curva- 
ture at the shock position, or possibly the shock curva- 
ture is infinite at the surface. A discontinuity in sur- 
face, or more correctly, bounding streamline curvature, 
would be expected to lead to a local pressure fall in a 
manner analogous to the incompressible flow over a 
curvature singularity. In the inviscid-flow problem 
postulated by Emmons, the pressure changes were 
found to occur in a very small region and led to solutions 
in which the shock pressure rise appeared to be less 
than that given by classical theory. From the analogy 
with the incompressible flow past a boundary curvature 
singularity, the streamwise extent of the pressure 
changes would be expected to depend on the strength 
of the singularity, which is dependent on the Mach 
number ahead of the shock and the local surface curva- 
ture. In a viscous flow, the phenomenon will clearly 
be associated with a local thickening or separation of 
the boundary layer, and is further complicated by the 
influence of the rapid pressure change itself on the 
equilibrium of the boundary layer. 

Downstream of a shock wave, the pressure changes 
must satisfy the boundary conditions imposed by the 
geometry of the surface, boundary layer, and wake, 
and, ia particular, the flow must return to the initial 
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free-stream conditions downstream. It is well known 
that for wholly subsonic flows, the pressure distribution 
over an airfoil at different free-stream Mach numbers 
can be related by simple compressibility formulas. <A 
significant development of this type of Mach number 
variation is its continuation in part to flows which in- 
clude a region of supersonic flow terminating in a shock 
wave. It is shown in reference 12 that the pressure dis- 
tribution some distance downstream of a shock can 
thus be approximately estimated from that in a wholly 
subsonic flow in the same way. 

From the above discussion of the transonic flow 
about an airfoil surface, it appears that the analysis of 
measured pressure distributions might be usefully con- 
sidered in three parts: (1) the supersonic flow upstream 
of the shock wave, (2) the shock pressure rise, and (3) 
the subsonic flow downstream of the shock wave. This 
approach is adopted in Section (2). 


(2) Analysis of Measured Pressure Distributions 


(2.1) Supersonic Region 


Theoretical and experimental studies of transonic 
flow have shown that the rate of change with free- 
stream Mach number of the local supersonic pressure 
distribution over an airfoil surface is small for free- 
stream Mach numbers near one. The typical family 
of measured pressure distributions of Fig. 1 illustrates 
that the supersonic pressure distributions are, in con- 
sequence, similar in form to the corresponding sonic- 
range (or J) = 1.0) pressure distribution. Analysis of 
this part of transonic pressure distributions has, there- 
fore, been based on the study of the differences between 
local pressure at high-subsonic values of 1) and that 
at the same chordwise position at 1) = 1.0. 

A prime requirement of the present scheme for the 
prediction of pressure distributions is the estimation 
of the pressure immediately upstream of a shock wave, 
pi. As suggested above, it is supposed that the major 
flow and geometrical parameters which determine this 
are included in the sonic-range distribution, so that the 
analysis can be restricted to the pressure (ratio) dif- 
ference (pi — Psr)/Ho. Values of this quantity ob- 
tained from a wide range of wind-tunnel tests are 
shown plotted against My in Fig. 2; on this, and the 
following Figures, each symbol denotes values for one 
airfoil configuration, and is defined in Table 5. Apart 
from the expected decrease in the pressure difference 
as Mo approaches unity, little consistency can be de- 
tected in this plot. It should be mentioned that no 
significance could be attached to shock position, except 
insofar as it is related to Mp for each airfoil configura- 
tion. However, one source of the scatter has become 
apparent from the results of tests on one airfoil section 
in all the high-speed tunnel configurations of the Aero- 
dynamics Division, NPL. It is found that the “‘solid- 
ity’’ of the Mach number freeze is greatly affected by 
the tunnel slotted-wall open/total area ratio used in the 
tests; moreover, the significance in this respect of a 


given open area is dependent on model incidence. To 
avoid confusion caused by these tunnel-wall effects, 
the analysis of the supersonic region is confined to re- 
sults obtained from the NPL 20-in. X 8-in. High-Speed 
Tunnel with 1/30th open/total area walls, and for 
zero incidence tests only. These test conditions were 
chosen for detailed study because they are known to 
give smaller interference effects in the transonic range 
than the other tunnel configurations used (see references 
13 and 14). In Fig. 3, these pressure differences are 
plotted against the quantity (WM) — M*)/(1 — M*), 
where J/* is defined as the free-stream Mach number 
for which local sonic speed is first reached at the air- 
foil surface crest,t on the assumption of shock-free 
flow. This parameter is introduced to represent dif- 
ferences in airfoil thickness and profile, and to assist in 
the comparison of results measured over different ranges 
of AJ. A useful mean curve can be drawn through 
the points of Fig. 3, and this is used to determine ~,/ [1 
for a given free-stream Mach number from the corre- 
sponding sonic-range distribution. 

In the following shock pressure-rise analysis, which 
is based on results obtained from all tunnel configura- 
tions, it is shown that p,//J) is not an important 
parameter in the determination of shock position. In 
the application of the present scheme for the calcula- 
tion of transonic pressure distributions, a value of p;//Zo 
is essential in the estimation of the supersonic region, 
and the use of Fig. 3 is suggested to give consistent, 
and approximately tunnel-interference free, results. 

To assist in the estimation of the supersonic pressure 
distribution upstream of a shock wave, an analysis, 
similar to that described above for p;/Ho, was carried 
out for the pressure differences at airfoil surface crests. 
Again, only results obtained in 1/30th open area tunnel 
walls with zero model incidence were used. Fig. 4 
shows these pressure differences, (Per Psr)/ Ho, 
plotted against (Jy) — /*)/(1 — M*); again, a useful 
mean curve can be drawn through the results. 


(2.2) Shock Pressure Rise 


It has been noted that discrepancies between meas- 
ured shock pressure rise ratios and the classical normal 
shock equation is one of the most disturbing features in 
attempts to predict transonic pressure distributions. 
Although Emmons has suggested a physical mecha- 
nism which can account for these differences, the addi- 
tional direct effects of shock-wave boundary-layer 
interaction appear to frustrate any attempt to relate 
this mechanism to observed flows. Shock pressure 
rises obtained from wind-tunnel tests are difficult to 
define accurately, and appear sensitive to the test 
conditions. Therefore, before an analysis of experi- 
mental results could be undertaken, it was essential 
to choose a consistent method of interpreting the data. 
The previously mentioned relation between the pres- 
sure distribution over an airfoil surface downstream of 


{ The position where the surface tangent is in the free-stream 
direction. 
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| | TABLE 1 
Sonic Range Pressure Distribution—Calculated by the Method 
of Reference 15 


x/e | 0.22 | 0.30 0.40 0.50 0.60 “0.80 0.90 1.00 
p/Ho 0.383 | 0.364 0.345 0.334 0.332 0.330 0.330 0.337 0.342 


\ TABLE 2 
~. Subsonic (Downstream of Shock) Part of Pressure Distribution, 
| > | Calculated From Measured Distribution for MJ, = 0.70 by the 
| i Glauert Compressibility Correction Factor 
| 0.70 0 75 0. 0 82 
P, 0.22 0.565 0.506 
0.30 0.588 0.530 
Empirical shock relation 0.40 0.618 0.568 0.508 
( — ———_—s— Exact normal! shock relation 0.50 0.645 0.600 0.550 0.527 
Curves for P, /H,=0'51,0-52 and 0-53 0.60 0.667 0.625 0.580 0 561 
0.90 0.730 0.699 0.669 0.657 
Fic. 6. Comparison of empirical, exact, and constant p»/H, 1.00 0.756 0.730 0.705 0.695 
shock relations. — 
TABLE 3 


Estimation of Shock Position 
(For the NACA 0009;, a = 2°, M* = 0.730) 


(a) Use of Figs. 3 and 4 


M M* — M*)/(1 — M*) — psr)/Ho (per — psr)/Ho ber/ Hy 
“3 0.75 0.02 0.074 0.075 0.054 0.437 
0.80 0.07 0.269 0.037 0.023 0.406 
0.82 0.09 0.334 0.029 0.017 0.400 
(b) Calculation of p;/Ho and pog/ Hy Loci 
————M, = = ———M, = 0.82———_—_— 
x/c pi/ Ho Pri Ho bi/ Ho pr / Ho Pr Ho 
0.22 0.458 1.115 0.511 
0.30 0.4389 1.183 0.519 
0.40 0.382 1.387 0.530 
0.50 0.371 1.427 0.530 0.363 1.455 0.529 
0.60 0.361 1. 462 0.528 
| PIH, distribution by interpolation 


from Peg)/Ho and 


| 
| 
locus from Pep|Ho and 
| 


(R-Poq)/Ho | | | 
from Pee /H, and Fig. 4 
| y Intersection gives shock 
| | 3 position | | 
| | | | 
P,/H, locustrom p,/H, 
Hp locus and Fig. 5. Sho | 
| | bution 
| 
| | 
| | 
| | 


Fic. 7. Illustration of scheme for predicting transonic pressure distributions. 
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a shock wave to the corresponding low-speed flow sug- 
gests a possible approach to this problem. It is shown 
in reference 12 that the Prandtl-Glauert rule gives a 
good representation of this type of flow behavior near 
the trailing edge, and this simple rule is used to define 
the pressure immediately downstream of observed 
shock waves. That is, in the analysis of wind-tunnel 
results, the pressure immediately downstream of the 
shock waves is defined as the value predicted at that 
position by the Glauert compressibility formula (applied 
to results for some lower Mach number) for the Mach 
number corresponding to the observed shock position; 
this pressure is denoted fog. This convention is re- 
garded as a means of interpolating in the neighborhood 
of the shock-wave in a consistent manner. To ensure 
consistency, only results which follow the Glauert 
compressibility variation at least at the trailing edge 
are used in the analysis. As would be expected, in all 
but a few of the examples considered po, differs from 
the measured, or more correctly, estimated values of 
ps, and it has not been possible to relate the differences 
to any flow or airfoil parameter. 

In Fig. 5 the shock pressure rise ratio pog/p; is shown 
plotted against p,//Zo, the results being obtained from 
a number of observed shock positions on sixteen dif- 
ferent airfoil configurations. It is seen that a reason- 
able correlation is obtained, independent of airfoil 
section or other flow characteristics. A mean curve is 
drawn through the points and the shock pressure rise 
relation given by this is used in the prediction of 
transonic pressure distributions described in the fol- 
lowing Section. It is interesting to note that this 
empirical relation is not greatly different from that 
obtained on the assumption that po,//7) has a constant 
value. This similarity explains the previously men- 
tioned insensitivity to p,/Hy of the comparison of 
values of px;/pi. An assessment of the significance of 
the scatter of the values on Fig. 5 in terms of airfoil 
shock position is given later. To assist in the later dis- 
cussion of the result, the mean curve of Fig. 5 is shown 
in Fig. 6 together with curves obtained with constant 
pe¢/ Hy and from the exact normal shock equation. 


(3) Application of the Empirical Scheme for the 
Prediction of Transonic Pressure Distributions 


In this Section, it is shown how the results of the 
analysis of Section (2) can be used to predict transonic 
pressure distributions from a knowledge of the geometry 
of an airfoil and its pressure distribution in wholly sub- 


Calculation of Supersonic Pressures Between Crest and Shock 


Moy = 0.80 Mo = 0.82 

x/¢ X  (p-psr)/Ho p/Ho x/c X (p-psr)/Ho p/H 
0.22 0 0.023 0.406 0.22 0 0.017 0.400 
0.30 0.34 0.028 0.392 0.30 0.28 0.020 0.384 
0.40 0.77 0.034 0.379 0.40 0.63 0.025 0.370 
0.455 1.00 0.037 0.375 0.50 0.98 0.029 0.363 


0.505 1.00 0.029 0.362 
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TABLE 5 
Summary of Experimental Results* 
Wind Model Inci- 
Tunnel Chord dence Sym- 
Airfoil (inches) (inches) Reference (degrees) bol 
NPL 491 36 X 14 9 ARC R&M 0 x 
3100 1/2 t 
1 + 
A 
2 > 
10°, RAE 36 X 14 5 Unpublished 0 A 
104 2 Vv 
6° RAE 20x 8 5 Unpublished 0 > 
104 2 =< 
4% RAE 20 X 8 5 Unpublished 2 V 
104 
0 < 
6% RAE 20X*8 5 Unpublished 2 ae 
102 4 
4% Bicon- 36 X 14 9 ARC 20,654 0 A 
vex 
NACA 20x 8 5 Unpublished 0 O 
0009; 2 A 
NACA 20 <X 8 6 Unpublished 0 € 
0009; 2 
with 
droop I 
NACA 20 X 8 6 Unpublished 
0009; 2 
with 4 
droop II 
9% 20 x8 5 Vickers- 2 
cambered Armstrong 
section (Aircraft) 
Ltd. 


Aero/FM / 
Report /002 


* All tests made in NPL high-speed tunnels. 


sonic flow. The method is described by proceeding 
through the stages of a typical calculation. 

The calculation of the upper surface pressure dis- 
tribution over a 9'/, per cent thick airfoil (NACA 
0009;) at 2° incidence for free-stream Mach numbers 
of 0.75, 0.80, and 0.82 is set out in Tables 1—4, and is 
illustrated for J) = 0.80 in Fig. 7. Table 1 gives the 
My = 1.0 (sonic-range) pressure distribution from 
crest to trailing edge as calculated by the theory of 
reference 15; this is also shown on Fig. 7. Table 2 
lists the pressure distribution measured at J/) = 0.70 
for the NACA 0009; at 2° incidence followed by the 
distributions calculated from this by the Glauert rule 
for My = 0.75, 0.80, and 0.82. The results for Mo = 
0.80 are also shown on Fig. 7. The next stages in the 
calculation are based on the empirical relationships 
established in the preceding Section of this paper. 
Table 3(a) shows the quantity (\J, — M*)/(1 — M*) 
tabulated against 1/4); the corresponding values of 
(pi — Psr)/Ho and (per — Psr)/H» are then read off 
from the mean curves of Figs. 3 and 4. In Table 3(b) 
values of p:/H» for each Mach number derived from 
Tables 1 and 3(a) are listed for the appropriate chord- 
wise positions. It should be emphasized that, on the 
basis of the preceding analysis, the difference (pi: — 
Psr)/Ho has significance only at the relevant shock 
position. However, in the construction of a solution 
for a given Mo, it will be seen to be convenient to de- 
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Present theory 


Experiment 
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Fic. 8. Comparison of results obtained in Tables 1 through 4 with experiment. 


rive a locus of possible values of p;///) for the pre- 
scribed Mo for a small chordwise region. The calcula- 
tion can be limited to a few points from the knowledge 
that the shock pressure rise relation leads to only small 
changes in with changes in To return 
to Table 3(b), the shock pressure rise ratio pog/p corre- 
sponding to each value of p;/H) is taken from Fig. 5, 
and px/Hpo calculated. In this way, a p2¢/Hp locus for 
each free-stream Mach number is derived from the 
oF corresponding ~;, //) locus. The broken lines on Fig. 7 
ey define these loci for .\/) = 0.80. Now, from the assump- 
tions in the analysis leading to the empirical relation- 
ships, the intersection of the pog/Hy locus with the 
downstream distribution estimated by the Glauert rule 
=, for the appropriate free-stream Mach number defines 
"3 the shock position for that Mach number; also at this 
; position, «s/c, the assumption in the analysis for p:/ Fo 
holds, and hence /; //) is determined. There remains 
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NPL 491 aerofoil 


NPL 431 aeroforl 
Fic. 9 (top). Comparison of theory expounded by Spreiter 


and Alksne? and experiment. Fic. 10 (bottom). Comparison 
of modified theory and experiment. 


the calculation of the pressure distribution between 
the crest and shock position. The pressure at the crest 
for each free-stream Mach number is obtained from 
the values of (ber — Psr)/Ho given in Table 3(a), and 
the sonic range value. To estimate the pressure dis- 
tribution from crest to shock position, that is, from 
Por/Hy to pi/Hp, it is sufficient to interpolate linearly 
with respect to chordwise position for the pressure 
differences from the sonic-range distribution. If 


A= (x Xer)/(Xs = Xcr) 
then 


(Pp — Psr)/Ho = (Per — Psr)/Ho + 
— Psr)/Ho — (Per — Psr)/Ho}-X 


In the example considered, this calculation is trivial for 
My = 0.75 and is only set out for My) = 0.80 and 0.82 
in Table 4. 

The estimated pressure distributions for My = 0.75, 
0.80, and 0.82 are shown on Fig. 8 together with the 
corresponding measured values. 


(4) Analysis of a Theoretical Solution 


To supplement the empirical analysis, some calcula- 
tions of transonic flows by the method of Spreiter and 
Alksne’ have been made. This theory, which leads 
to an iterative solution of an integral form of the 
transonic flow equations, is the most satisfactory 
theoretical method available. Fig. 9 shows the pres- 
sure distributions calculated for a 4 per cent thick 
airfoil (NPL 491) at zero incidence compared with those 
obtained by experiment. It is seen that while the 
theoretical results exhibit many of the features of the 
transonic flow about an airfoil surface, the predicted 
shock positions are significantly downstream of those 
observed at the same free-stream Mach number. 

It is suggested that the large discrepancies in shock 
position can be attributed to the assumption in the 
theory that the pressure rise in the vicinity of the shock 
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is that given by the normal shock equation. It has 
been suggested earlier that the presence of a normal 
shock on a curved surface leads to a reversal of the 
pressure gradient normal to the surface, and to a dis- 
continuity in bounding streamline curvature and a 
consequent pressure fall. Now, in the application of 
the theory, no solution is obtained at the shock position 
itself, and also, the boundary condition assumed on the 
airfoil surface is incompatible with a reversal in the nor- 
mal pressure gradient. Thus, the theory permits no rep- 
resentation of local flow behavior at the shock position. 
The empirical shock relation established in Section (2) 
would, therefore, seem to be more consistent with the ap- 
proximations made in reference 7, as this is defined to in- 
clude the overall pressure changes due to shock-boundary 
interaction. The calculated solutions can easily be 
modified to include the empirical shock relation as 
follows. 

First, it has been found that the calculated pressure 
distributions downstream of the shock positions are 
similar to those estimated by the Prandtl-Glauert rule 
applied to a wholly subsonic solution. This is in ac- 
cordance with the previously noted behavior of meas- 
ured transonic distributions. Therefore, it seems 
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Fic. 11. Comparisons of predicted and measured shock-wave 
positions. 
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Fic. 12. Effects of leading-edge geometry modifications on 


pressure distributions. (a) NACA 0009; and modifications, 
My) = 0.75, a = O°. (b) Goldstein 14/47 and modifications, 
M, = 0.50, a = 5°. 


reasonable to extend the subsonic part of the calculated 
pressure distributions forward of the original shock 
positions by use of the Prandtl-Glauert rule. Second, 
the physical mechanism of the development of the 
supersonic part of the flow at a prescribed free-stream 
Mach number suggests that the associated surface 
pressure distribution is little affected by shock position, 
and the analysis of Section (2.1) confirms this view. 
The supersonic distribution calculated by the method 
of reference 7 may thus be taken as a locus of possible 
values of the shock upstream pressure. Then, in a 
similar way to the construction of transonic pressure 
distributions described in Section (3), the empirical 
shock relation defines a fog locus, whose intersection 
with the appropriate subsonic distribution determines 
the shock position. This procedure has been applied 
to the results shown in Fig. 9 to give those shown in Fig. 
10; it isseen that good agreement with experiment is then 
obtained. To illustrate further the comparison of both 
the original and modified solutions with experiment, the 
shock positions for the NPL 491 airfoil, and two other 
examples, are shown in Fig. 11 as functions of Wo. 


(5) Discussion 


In this Section, the significance of the analysis of 
Section (2) is considered in relation to the physical 
mechanism of transonic flows. The applicability of 
the method for calculating transonic pressure distribu- 
tions is discussed, and the accuracy of the results ob- 
tained is assessed by comparison with experimental 
data. 
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Sketch A 


(5.1) Physical Mechanism of Mixed Subsonic /Supersonic 
Flows 

Before discussing the significance of the empirical 
relations established in this paper, the following con- 
vention will be adopted. It is noted in Section (1) that 
the pressure fall associated with the normal pressure 
gradient reversal through a normal shock wave on an 
airfoil surface is distributed over a limited chordwise 
region, and that this effect is probably influenced by 
the shock-wave boundary-layer interaction in a real 
flow. As it is not possible to calculate this pressure 
distribution reliably, it is convenient to suppose that 
the total pressure change can be represented as occur- 
ring at the shock position. This convention is illus- 
trated in Sketch A. The full line represents the true 
pressure distribution and the hatched area the pres- 
sure changes due to the curvature singularity following 
the shock, or to some other aspect of shock/surface 
interaction. The assumption that this pressure dis- 
tribution is concentrated at the shock position is 
equivalent to the replacement of the true distribution 
downstream of the shock by that shown by the broken 
line. 

It is suggested that the empirical relations estab- 
lished are consistent with the attainment of the equiva- 
lent shock-free flow* downstream of a shock wave. 
The pressure distribution downstream of a shock given 
by the convention defined above would thus be that 
obtained on the assumption of shock-free flow. This 
hypothesis accounts for the finding that measured 
pressure distributions in the vicinity of the trailing 
edge of an airfoil are related by a subsonic flow com- 
pressibility correction formula, even though shock 
waves are present in the flow upstream. Also, the 
relation for the shock downstream pressure defined on 
the same basis (as po¢/ fo) established in Section (2.2), 
suggests that the equivalent shock-free flow has some 
significance at this position. It is suggested that for 
a given shock upstream pressure (/;///o), the difference 
between the value of p2/Hp given by the exact normal 


* The “equivalent shock-free’” flow is that obtained for the 
prescribed free-stream Mach number on the assumption of sub- 
sonic flow, although it will usually include a (unreal) region of 
supersonic flow; here it is taken as that given by the Glauert rule 
applied to a wholly subsonic flow. 
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shock equation and that estimated on the assumption 
of shock-free flow (p2¢/ H/o) is the total effect of the pres- 
sure changes due to the shock/surface interaction. 
The physical mechanism of the flow can then be de- 
duced as follows. The pressure changes through a 
shock wave are to be thought of as including a full 
normal shock pressure rise followed by a pressure fall 
dependent on the shock/surface interaction. The 
shock position for a given free-stream Mach number 
is then determined by the necessity for the resultant 
shock downstream pressure to be continuous with the 
shock-free flow consistent with the prescribed free- 
stream conditions. This hypothesis implies that the 
replacement of an artificial, smoothly decelerating, 
supersonic flow by a real supersonic flow terminating 
in a shock wave does not significantly affect the sub- 
sonic flow downstream. 

Attention is drawn to the measured pressure distri- 
butions shown in Figs. 12(a) and 12(b). In Fig. 
12(a), wholly subsonic flow results for the basic NACA 
0009; section are compared with those at the same 
Mach number and incidence for two modified forms 
of this section. The modifications consist of different 
forms of dropped leading-edge extension which give 
large changes from the basic pressure distribution 
forward of the crest of the surfaces, but are seen to 
have a negligible effect over the rear part of the sec- 
tions. Fig. 12(b) shows a similar comparison between 
the subsonic-flow pressure distributions over a 14 per 
cent thick airfoil with two leading-edge shapes; again, 
although there are large differences in the pressure dis- 
tribution near the leading edge, there is no effect of 
the change in leading-edge shape on the rest of the 
pressure distribution. The significance of the change 
in the type of supersonic flow over the forward part of 
an airfoil mentioned above is thus suggested as analo- 
gous to a change in airfoil shape under conditions in 
which the flow type is unchanged, in that neither have 
a significant effect on the downstream flow. In either 
case, it would seem necessary that the flow or geometry 
change has no more than a local effect on the boundary- 
layer development, particularly for flows with circula- 
tion. 

The hypothesis that the shock pressure rise ratio 
pi represents the overall pressure changes associ- 
ated with a transonic shock wave would seem to be 
supported by the close agreement with the empirical 
curve for of Emmons’ theoretical results, as 
shown on Fig. 6. However, the strength of the 
boundary curvature singularity associated with a nor- 
mal shock wave is simply demonstrated by Emmons to 
be related to the shock upstream pressure and to the 
local surface curvature. Hence, with each possible 
shock position the overall shock pressure changes 
should be related to both the local p,//7/) and surface 
curvature; whereas the present analysis for the overall 
shock pressure ratio gives a result which is independent 
of surface curvature. It thus appears that the dis- 
crepancy between exact and observed shock pressure 
rises cannot be entirely explained by consideration of 
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Fic. 13. Comparison of 


expb. 


calculated and observed shock 


positions. (Ringed values were not used in analysis for shock 


pressure rise. ) 


an effective surface curvature singularity, at least for 
viscous flows. 

Another interesting feature of the shock relation is 
illustrated in Fig. 6; it is seen that, throughout the range 
of Hp considered, Hp varies little, and that a good 
approximation would be to assume that pog = 0.528 
HH, which corresponds to sonic conditions. It thus 
appears that shock waves occur at approximately the 
downstream sonic point of the corresponding shock- 
free flow. This would be in accordance with the well- 
known suggestion that a disturbance to a closed shock- 
free supersonic flow region propagates in such a way as 
to lead to a breakdown of the smooth supersonic to 
subsonic deceleration at the downstream end of the 
supersonic region. Of course, the establishment of a 
shock in this way must lead to change in the supersonic 
flow field, and it cannot be argued’ directly that the 
shock position must coincide with the corresponding 
shock-free flow sonic point. However, such a flow 
mechanism would accord with the results of the present 
analysis. 


(5.2) Applicability of Empirical Scheme 


In the analysis leading to the empirical relationships 
described in Section (2) only selected experimental re- 
sults are used. Measured pressure distributions which 
appear to be incompatible with the physical basis of 


the analysis are excluded. First, only results for which 
the conception of a Mach number freeze is approxi- 
mately valid have been used, so that pressure distri- 
butions which include leading-edge suction peaks are 
not included in the analysis. Such flows are closely 
related to airfoil leading-edge geometry, and cannot 
be simply related to the sonic range distribution by a 
function of free-stream Mach number only. Also, in 
the analysis of supersonic regions, empirical relations 
are derived from results obtained under one set of wind- 
tunnel conditions. This restriction was necessary to 
avoid the complicated and not fully understood effects 
of slotted-wall tunnel interference. The second major 
restriction follows from the definition of the shock 
downstream pressure in terms of the equivalent shock- 
free flow. To ensure consistency in the analysis based 
on observed shock positions, only results which could 
be shown to be coincident with the corresponding shock- 
free flow at least at the trailing edge have been con- 
sidered. The analysis has also been limited to pres- 
sure distributions between the crest and trailing edge 
of airfoil surfaces, primarily because it is only possible 
to estimate the sonic range distribution for round- 
nosed airfoils at incidence over this region. Also, no 
separate consideration of the lower surface of an airfoil 
at incidence is included. In cases when the flow is of 
mixed subsonic/supersonic type, the present analysis 
is relevant; if the lower surface flow is wholly subsonic, 
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the pressure distribution can again be estimated ap- 
proximately by the Glauert rule. 

To assess the value of the present scheme for esti- 
mating transonic pressure distributions, it has been 
applied to a large number of cases for which experi- 
mental results are available. The comparison of theory 
and experiment is shown in Fig. 13 in terms of shock 
position; many examples not used in the analysis are 
included. In nearly all cases, good agreement is ob- 
tained, even though in some of the examples the above 
restrictions were not satisfied. 

The present analysis is based on shock-free pressure 
distributions obtained from measured low-speed dis- 
tributions, which ensures that major viscous effects 
are included throughout. However, recent applica- 
tions'® based on calculated inviscid flow shock-free 
pressure distributions, have shown that useful results 
can be obtained without further reference to experi- 
mental results. 
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these values of m, moment and normal force diagrams 
have been determined, and are shown on Fig. 2. The 
diagrams for linear elastic (or linear viscous) material 
(case m = 1) are also shown. As expected, the bend- 
ing moments (especially the peaks at the end) are con- 
siderably reduced as a result of the nonlinearity. On 
the other hand, the length over which bending mo- 
ments are significant is considerably increased. 

To check the assumed shape of the deformed shell, 
some preliminary tests have been performed. The 
creep deformation of a lead cylinder with clamped end 
(2 in. internal diameter, 0.2 in. wall thickness) coincided 
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fairly well with that used in the theory (for m ~ 3), as 
indicated in Fig. 3. 
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Initially Warped Sandwich Panel Under 
Combined Loadings' 


CHIEH C. CHANG* ann BERTRAND T. FANG** 


Unwersity of Minnesota 


Summary 


Based on small deflection theory, differential equations for the 
elastic bending of an orthotropic weak-core sandwich panel with 
small initial warping are derived by the variational energy 
method. The applied loads consist of arbitrarily distributed 
transverse loads and eccentrically applied edge loads and/or 
edge moments. For the case of a simply supported rectangular 
panel, solutions of the differential equations are obtained in the 
form of double Fourier series. As a practical example, numerical 
values of the maximum stresses and deflections for a square sand- 
wich panel are calculated and presented in the form of graphs. 

A simple design criterion is suggested for considering approxi- 
mately the effect of initial warping and transverse pressure. 
The method is particularly useful for preliminary design. 


Symbols 

A = surface area of panel 

B; = (1+ + ratio of fac- 
ings’ bending rigidity to that of sandwich panel 

D = bending rigidity of plate 

E = modulus of elasticity of facing 

= strain tensor 

g = 

Ges = moduli of rigidity of core 

Bus = lateral dimension of panel 

m = 

Mag = intensity of external moments applied at panel 
edge 

Nag = loading intensity applied at panel edge 

Po = wE't'l?/(L,~1 + m)(1 — yw*)], one quarter of 
critical cylindrical load (L,— = ) 

q = intensity of transverse loading acting on panel 
surface 

Fs = P,/iG,., core shear stiffness parameter for sand- 
wich panel 

s = boundary of panel edge 

t = thickness of facing 

i = thickness of core 

» = displacement in the x; direction 

U,V, wi 

u°s = initial deflection of panel 

us! = u°3 + u;, total deflection of panel 

U = potential energy of the system 

= rectangular Cartesian coordinates 

B = L,/L,, length-width ratio of panel 

6 = variational sign 

65; = Kronecker delta 

Ne = components of outward unit normal vector to 


boundary 
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m = Poisson's ratio of facing 
Tij = stress tensor 
Subscripts: 
i,j,k, 1 take on the values 1, 2, and 3 
a, B, y,6 take on the values 1 and 2 
Superscripts: 


convex facing 
concave facing 


” 


core 


Introduction 


__—- PANEL is composed essentially of a rela- 
tively thick core of low density which is bound 
between two thin sheets of strong material. To an- 
alyze the elastic behavior of such panels the usual thin 
plate theory cannot be applied because (a) the sand- 
wich panel is a composite structure where the homo- 
geneity assumption does not hold, (b) shearing defor- 
mation produces appreciable effect in a sandwich panel 
and can no longer be neglected, and (c) the elastic 
property of the sandwich core is usually orthotropic. 
A theory has been developed by Libove and Batdorf! 
to describe the approximate elastic behavior of a sand- 
wich panel by the so-called thick plate theory using 
effective elastic constants and taking into considera- 
tion the shearing deformation and the orthotropy of 
the panel. Other investigators*~‘ consider separately 
the different behavior of the different layers of the sand- 
wich panel and take into consideration the core shear- 
ing deformation. They assume that in the core plane 
sections initially perpendicular to the panel remain 
plane under deformation, but rotate by an amount pro- 
portional to the slope of the transverse deflection of the 
middle surface. It was shown’ that this procedure 
would not yield an exact result unless the panel was 
simply-supported at the boundary. Hoff** used a 
variational method to analyze certain simple sandwich 
structures of the so-called weak core. His method has 
the advantage of simplicity and elegance in obtaining 
the governing differential equations and boundary 
conditions. Also, his result is more reliable than those 
of other works cited above because he has only assumed 
that in the sandwich core plane sections originally per- 
pendicular to the panel remained plane during deforma- 
tion, and he has made no a priori assumption as to the 
amount of rotation. The extension of Hoff’s method 
to more complicated sandwich constructions was done 
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Fic. 1. Initially warped sandwich panel under specified loadings. 


recently by Chang and Ebcioglu’ and Chang, Fang, 
and Ebcioglu.® 

The present paper is an extension of the work of refer- 
ence 9 on the theoretical investigation of the elastic 
bending and buckling of sandwich panels. The sand- 
wich panel under consideration is rectangular in shape 
and is initially warped. The amount of warping is 
restricted to within the framework of small deflection 
theory. The initial warping is assumed given and may 
be due either by initial imperfections or to a temperature 
gradient across the panel, or to both. The two facings 
of the panel are assumed isotropic but may be made of 
different materials and may have different thicknesses. 
The core of the panel is elastically orthotropic, and is 
of much lower stiffness than the facings. The panel is 
subjected to an arbitrarily distributed transverse load 
acting on the concave side of the panel, and arbitrarily 
distributed forces and moments along the edges of the 
panel. Fig. 1 shows a rectangular sandwich panel 
under specified loadings and the orientation of coordi- 
nate axes with respect to the panel. 


Assumptions 


The simplifying assumptions made in the analysis are 
briefly outlined in the following. (For a more complete 
discussion of the assumptions, see reference 8.) 

(a) The forces in the plane of the panel are taken up 
entirely by two facings of the sandwich panel. The 
facings obey the classical thin plate theory, with the 
stretching of their middle planes being of particular 
importance. 

(b) The core behaves like a homogeneous orthotropic 
elastic continuum. The three planes of elastic sym- 
metry are the coordinate planes parallel to the surface 
and the edges of the panel. The core carries the main 
part of transverse shear loading and undergoes con- 
siderable shear deformation. As a consequence of 
assumption (a), this transverse shear is constant 
through the thickness of the core. 

(c) The deflections of the panel are small, and are 
constant through the panel thickness. This precludes 
symmetrical wrinkling of the facings. 

(d) Failure of the panel does not occur at the bonds. 

Assumptions (a) and (b) are based on the relative 
stiffness of the core and the facings. They are on the 
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conservative side for the so-called weak-core* sandwich 
structure. Assumption (c) imposes a restriction on the 
magnitude of the transverse loads allowed. Symmet- 
rical wrinkling of sandwich facings is not likely to 
occur in a panel under combined bending and com- 
pression, and will not be treated here. Similar as- 
sumptions have been used previously by Leggett and 
Hopkins,” Hoff,’ Chang and Ebcioglu’ and others in 
the analysis of sandwich structures. The satisfactory 
agreement found between experimental data and theo- 
retical resulsts indicates that these assumptions are 
justified. 


Displacements Characterizing the Deformation 
of a Sandwich Panel 


Since the convex and concave facings of the sandwich 
panel obey the law of the usual plate theory, their be- 
havior can be completely described by the displace- 
ments w’(X), %2) and uw” (x1, of their middle planes. 
The displacements of the core are denoted by #;(x, 
X2, X3). According to assumption (c), the deflection of 
the panel is constant through the panel thickness, or 
= u"3 = = U3(X1, Furthermore, from as- 
sumption (b) that the transverse shear is constant in 
the core, and from the continuity of displacement at 
the interface of core and facing, the displacements of 
the core #,(X1, X2, ¥3) can be written approximately in 
terms of the displacements of the facings u’, and 7”, as 


ig = — u"Q)/U (x3 + =(a=1,2 


This equation implies the approximation that the dis- 
placements at the interface of the core and either of 
the facings are equal to the corresponding displacements 
at the middle plane of that facing. From the above 
discussion, it is seen that in order to describe the de- 
formation of the sandwich panel, it is sufficient to know 
the five displacements x2), and 
3(X1, X2). Of these displacements uw’, and u”, are not 
independent. Notice that the displacements «’, and 
wu", give rise to two pairs of equal and opposite normal 
and tangential resultant forces in the facings. These 
two pairs of forces form a bending and a twisting couple 
which oppose the bending and twisting of the sandwich 
panel as a whole. This condition gives a relationship 
between uv’, and w”,. In this analysis we shall assume 
that the Poisson’s ratios for the two facings are equal, 
which is approximately true for practical face mate- 
rials. Under this condition, the relationship between 
uw’, and u”, can be reduced to a simple form 


u", = = —mu', (1) 


where m is defined as E’t'/E”t".+ Therefore the de- 
formation of the sandwich panel can be completely 


* By weak-core we mean that the core is weak in carrying 
direct stress. It does not necessarily mean that the core is weak 
in carrying transverse shear. 

7 This relationship is introduced here primarily for the sake 
of convenience. It replaces the force equilibrium equations in 
the plane of the panel which would otherwise be obtained later. 
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described by the three displacement components u’, 
(or u”,) and u;. It should be noted that wu; is the de- 
flection measured from the initially warped configura- 
tion of the panel and the wu’, are the displacements 
associated with the bending of the panel. No coordi- 
nate is provided for the displacements due to the edge 
compression and shear. The relevant material con- 
stants which characterize the sandwich panel are 
Young’s moduli E’ and Poisson's ratio = uw” 
of the isotropic convex and concave facing, and the 
shear moduli G,, of the orthotropic core. The differ- 
ential equations which relate the three unknown dis- 
placements uw’, and u; to the given loads g and Ny 
the initial deflection «°;, and the known material con- 
stants EZ’, E”, uw and G,, will be derived in the following 
sections. The derivation is based on the principle of 
minimum potential energy. First the energy integrals 
for the deformed panel and load system are formulated. 
Then the variational methods are used to minimize the 
total potential energy expression. This minimization 
process gives the differential equations of the panel 
together with associated natural boundary conditions. 


Total Potential Energy of the Panel and Load 
System 


On the basis of our assumptions, the total potential 
energy of the sandwich panel and the load system was 
obtained in reference 9 as follows :* 


= + m) f BU, + x 
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As it stands, this expression is quite general: it does 
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* ie this and the equations following, repeated Greek or Latin subscripts, e.g., yy, will imply summation. 
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Resultant forces and moments acting on a sandwich 
panel. 


FIG 


not specify a shape for the panel and it does not restrict 
the way in which applied loads are distributed. 

In obtaining the above expressions, we have adopted 
the assumption inherent in the small deflection theory 
that the forces .V,, in the neutral surface of the panel 
remain unchanged during transverse deflection of the 
panel. This is reflected in the fact that we have not 
included in Eq. (2) the strain energy of the stretching 
of the neutral surface of the panel due to direct-edge 
compression and shear, and in the fact that it is per- 
missible to write the potential energy of the edge forces 
in terms of the interior forces and displacements in the 
form of the last integral in Eq. (2). 

The energy expression Eq. (2) will be minimized with 
respect to the virtual displacements 6u; to give the 
differential equations of the panel together with asso- 
ciated natural boundary conditions. This same energy 
expression can also be used to find solutions by approxi- 
mate methods such as the Rayleigh-Ritz method’ 
for the upper bound and the Lagrangian multiplier 
method'! for the lower bound of the solution. 


Differential Equations and Associated Natural 
Boundary Conditions 


For a state of stable equilibrium, the total potential 
energy of the system assumes a minimum,'* and from 
the principle of the variational method it is necessary 
that the first variation of the potential energy shall 
vanish for any virtual displacements 6z;. 

The first variation of the potential energy can be 
shown to be*® 
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In order that the first variation of the total potential 
energy vanish identically, every integral in Eq. (3) 
must be individually zero. Since the virtual displace- 
ments 6u, in the surface integrals are entirely arbitrary, 
the expression in each individual brace under the sur- 
face integral sign must vanish identically. Therefore 


— pty — + — 
- m 
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Note that in (4) there is no summation over a. These 
are the three differential equations of equilibrium of 
the sandwich panel. It is desirable to examine these 


782 JOURNAL. OF THE AEROSPACE 


SCIENCES—OCTOBER, 1960 


equations from equilibrium considerations. A differ- 
ential element of length dx;, width dx. and height 7 + 
t’ + t” is cut out from the panel. Figure 2 shows such 
an element and the resultant forces and moments act- 
ing on it. Considered as a rigid body, this element 
has six degrees of freedom and is governed by six equa- 
tions of equilibrium. Eq. (4) derived above repre- 
sents the two moment equilibrium equations in the 
Xq — X3 plane. The remaining moment equilibrium 
equation is satisfied identically by the symmetry of 
stress tensor in the facings. Eq. (5) is the force 
equilibrium equation in the transverse direction. The 
other two force equilibrium equations reduce to Eq. (1) 
obtained earlier as a relationship between the displace- 
ments in the concave facing and the convex facing. 


The vanishing of each individual line integral in Eq. 
(3) furnishes the natural boundary conditions of the 
panel. It is sufficient that, along the boundary edges 
of the panel, 


Et’ M 
2) + (1 — + + 7 9 or 6u’, = 0 
+ 
2 [(1 U3 + = 0 or =  } (6) 
D’ + D” [lt+m , 
> + (4 We ina) + + = 0 or bu; = 0 


The conditions on the left-hand side of Eq. (6) imply 
that the limiting values of the internal forces and mo- 
ments, as the edge of the panel is approached, must be 
in equilibrium with certain prescribed forces and mo- 
ments externally applied at the edge of the panel. 
Notice that these boundary conditions correspond to 
the three edge conditions in E. Reissner’s theory'® of 
plates with transverse shear deformation considered. 


The set of Eqs. (4) and (5) are three coupled linear 
partial differential equations in terms of the three 
unknowns and It should be kept in mind 
that the NV, appearing in Eq. (5) are the interior forces 
acting in the neutral surface of the panel. Only in the 
case where the externally applied edge compression 
(tension) and shear are uniformly distributed, do these 
interior forces N,g, remain the same throughout the 
panel and equal their external values. When the edge 
forces are not uniformly distributed, the interior forces 
Ng depend on their location in the panel and can be 
determined by finding the Airy’s stress function in the 
plane stress problem. In such cases, however, Eq. (5) 
above becomes a linear partial differential equation 
with variable coefficients, and the solution of Eqs. (4) 
and (5) may become quite difficult. 


Therefore, the behavior of a deformed sandwich panel 
is completely characterized by the three displacement 
components w’, and uw; which are solutions of the linear 


partial differential Eqs. (4) and (5) subjected to the 
boundary conditions Eqs. (6). When these displace- 
ment components are known, the stresses in the panel 
and the buckling load of the panel can easily be calcu- 
lated. The solution given in the next section for a 
simply supported sandwich panel will serve to illus- 
trate these statements. 


Solution for a Simply Supported Sandwich 
Panel 


(a) Local Displacements and Deflection 


The solution for the displacement components of a 
rectangular sandwich panel simply supported at the 
edges and subjected to arbitrarily distributed trans- 
verse load and uniformly distributed edge compression 
can be obtained by expanding the displacement com- 
ponents, the transverse pressure, and the initial de- 
flection in terms of the eigenfunctions of the region. 
The solution as given in reference 9 is rearranged, with 
the introduction of appropriate dimensionless param- 
eters, to appear as follows: 
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= > (Kym? + + = | COS an — 
L( 1 m) m=1 n=1 Bivins 0 Ly 
Vv + A (K zm* + Kyn w mn Ymn | Sin I cos 
LA + m) m=1 n=1 mn +z “y 
| (Kom: + mn sin —— 
=] 0 L, Ly 


In these equations w®,,,, and Gm, are, respectively, the coefficients of the Fourier expansion of the initial deflection 
w"(x, y) and the transverse pressure q(x, y), and 


Cmn = Ve + n's*) em? | +1 


gret(m® + n%8°)? + 1, | + n'a?) +2 +1 


Amn = + n*B’)? — Kzm* — Kyn*p*| + (m* + re(gm* + + 

The appropriate dimensionless parameters are defined as follows: 6 = L,/L, is the aspect ratio of the panel, and 
r, = Po/G,,lis the core shear stiffness parameter of a sandwich panel, roughly proportional to the ratio of the flexural 
rigidity to the shearing rigidity of a sandwich panel. The ratio g = G,./G,, is the ratio of the shearing rigidity of 
the core in the two directions. By, = 7*(D’ + D”)/PoL,’ is a measure of the importance of the bending rigidity of 
individual faces relative to that of the whole panel. It is usually a small quantity, and taking its value as zero is 
equivalent to saying that the facings of the sandwich panel behave as membranes. A, = N,,/Po and A, = N,,/Po 
are dimensionless edge-compressive forces, where Py = 7°E’t'??/L,°(1 + m)(1 — yp’) is one fourth of the critical buck- 
ling load per inch width in a sandwich panel of infinite width (L, ~ ~) with shear deflection and bending rigidity of 
the facings neglected. 


(b) Stresses in Facings and Core 

The stresses in the facings of a sandwich panel are composed of two parts: 

(1) The stresses due to direct edge compression. In each facing these stresses are uniform compression; they are 
not the same, however, in both facings, being proportional to the respective modulus of elasticity of each facing. 

(2) The stresses due to bending of the sandwich panel. These stresses can again be decomposed into two parts, 
i.e., 

(i) The stresses due to bending of the panel as a whole about the neutral surface. These stresses are represe ted 
by those of the middle planes of the facings and will be called membrane stresses hereafter. They can be calculated 
from the stress-strain relationship for the plane stress problem, with the strain components being those associated 
with the middle plane displacements u’, and u”,. The membrane stresses in the concave facing are: 


m=1 n=1 Amn 


Po = mn + n 


t"l m=1 n=1 Amn 


| Kam’ + K,n’B?)w mn + sin 


= , 2 2 
= mnB(c ma tC na) | + K,n'8? += = = | max 
m=1 L, L, 
The membrane stresses in the convex facing are —t”/t’ times those in the concave facing. 

(ii) The stresses due to the bending of each individual facing about its own middle plane. We shall call these 
stresses curvature stresses hereafter. These stresses vary linearly through the facing thickness and can be calcu- 
lated from the transverse deflection w according to the usual plate theory. The stresses at the top edge of the con- 
cave facing are 


PA1 + 8) (m* + ms L? mmx 
= om? (Kym? + mn + =P, sin sin 
(Kym* + mn + =P, | sin L, 
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The corresponding stresses in the convex facing are m times as large. 

The resultant stresses in the facing are the sum of the uniform compression, the membrane stresses, and the curva- 
ture stresses. The top surface of the concave facing is the most stressed part of the facings because in addition to 
the uniform compression, the normal components of the membrane and curvature stresses are also compressive 
stresses there. 

The stresses in the core of a weak-cored sandwich panel are predominantly shearing stresses. They are given by 


the following expressions: 


given above are good design criteria. 


Tx = > K mn + 2D | SM cos 

Tyz = + mn + 2D cos sin 

iL, m=1 n=1 Amn 0 

where C4mn = m(m? + n7B") gr,(m? + + 1 


1— 
Cnn = nB(m* + n°B*) E + + 


Because of our assumption about the mode of deformation of the core, these values are average shearing stresses and 
are somewhat smaller than the actual maximum shearing stresses. However, usually core cells are so thin that 
failure would occur as the result of shear buckling rather than yielding, and under such circumstances the stresses 


I (c) Determination of Critical Edge Compression 
It is seen that A,,,, appears in the denominator of the displacements and stresses. When A,,,. = 0, the stresses and 
deflections become infinite. This corresponds to the buckling of the panel. A series of buckling loads are obtained 
. from the solution of the equation A,,,, = Oas 
K,m? + K,n’p? = 
3 9 1 — 9 999 
+ | r.(gm® + n°B*) + 1 
4 + + re | | m? + + g + — — + 1 
To each pair of the integers m and there correspond several combinations of the buckling loads K, and K,. For 
_ any given value of the ratio A,/A, there is a pair of integral values of m and n which makes K, and K, a minimum. 
” This is the critical buckling of the panel that actually occurs and these values of m and n give the actual buckling 


pattern of the panel. 


Numerical Example of an Initially Warped 
Square Panel 


Numerical calculations were carried out for the maxi- 
mum stresses and deflections arising in a simply-sup- 
ported and initially warped square sandwich panel 

$ which is subjected to a uniformly distributed com- 
pression NV, on two opposite edges and a uniform 
transverse pressure g on its concave surface. The 
initial warping of the panel is represented by w® sin 
(rx/L,) sin (ry/L,), which is a satisfactory approxi- 
mation for certain important cases such as the warping 
of a panel due to a constant temperature gradient 
across the panel. 

Since the effects of the initial deflection of the panel 
and of the transverse pressure are superimposable 
within the scope of small deflection theory, the stresses 
and deflection caused by these two different effects 


could be calculated separately. Because the initial 
warping is approximated by the first term of its Fourier 
series expansion, the stresses and deflection due to the 
combined action of the edge compression and this initial 
warping are also given by a single term and therefore 
can readily be calculated. On the other hand, since 
the uniformly distributed transverse pressure is repre- 
sented as a double infinite series, so also are the stress 
and deflection. The numerical evaluation of these 
series may involve some difficulty because the location 
of the maximum stresses and deflection may be un- 
known and because these series may converge rather 
slowly. For the square panel considered, and under 
the specific loading condition, the maxima of the de- 
flection « and the normal stresses in the facings r,,;, 
Tyy, tr2, and #,, occur at the center of the panel. Their 
series representation alternates in sign with respect to 
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both summing indices, and the summation of a limited 
number of terms gives an adequate degree of accuracy. 
The shearing stresses in the core 7,, and 7,, vanish at 
the center of the panel and reach their maxima at the 
midpoints of the edges of the panel. But the series 
representation for these stresses alternates only with 
one summing index and converges more slowly. The 
computations were done on a Univae 1103 digital 
computer, and partial sums of the series for these 
maximum stresses and deflections were evaluated for 
the first 256 nonvanishing terms of each series. The 
facings of the panel are also subjected to shearing 
stresses 7,, and #,,. These stresses vanish at the center 
and are maximum at the corners of the panel. The 
series representation for these stresses does not alter- 
nate with either summing index, and converges slowly. 
Fortunately these shearing stresses in the facings are 
much smaller than the corresponding normal stresses. 
Furthermore, at the location where the normal stresses 
are maximum, the shearing stresses vanish. There- 
fore, these shearing stresses in the facings are not im- 
portant for design purposes and are not computed. In 
addition to the stresses due to the bending of the panel, 
there is, of course, uniform compression in the facings 
due to direct edge compression. The above statement 
about the location of the maximum stresses and de- 
flection may not apply for panels with unequal length 
and width or for panels under other loading conditions. 

The range of panel parameters employed in the cal- 
culation of the stresses and deflections were r, = 0, 
0.1, 0.2, 0.3; g = 0.8, 1.0, 1.25; » = 0.3 and B, = 0. 
As is known, the shear stiffness parameter 7, is a meas- 
ure of the relative importance of the shear rigidity and 
the bending rigidity of a sandwich panel. Taking 
r, = 0 amounts to saying that the panel is perfectly 
rigid in shear, and the sandwich panel is reduced to the 
usual thin plate with an equivalent bending rigidity. 
As r, increases, the sandwich panel behaves less and less 
like the usual thin plate, it weakens owing to the re- 
duced shear rigidity, its deflection increases, and its 
buckling load decreases. In the calculations B, was 
taken as zero, which seemed to mean that we were 
neglecting the flexural rigidity of the facings of the 
sandwich panel in comparison with the flexural rigidity 
of the whole panel. In consistence with this assump- 
tion, it would seem unnecessary to calculate the curva- 
ture stresses in the facings. That this is not so can be 
explained by the fact that the resisting moment of the 
sandwich facings is of a smaller order of magnitude 
than the maximum curvature stresses due to this 
moment. A similar situation is encountered in the 
theory of plates where we start by neglecting the 
transverse shear deformation but end up finding the 
transverse shear stress. 


Graphical Representation of Results 


An obvious way of presenting graphically the nu- 
merical results would be to plot nondimensionally the 
maximum stresses and deflections against the edge 
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compression with r, and g as parameters. Since the 
panel parameters r, and g range from 0, 0.1, 0.2, to 0.3 
and from 0.8, 1.0, to 1.25, respectively, there would be 
ten curves for each stress or deflection (for 7, = 0, one 
curve would represent all values of g). It is found, 
however, that constructing the graphs in a different 
way is more advantageous. We shall define a reduced 
edge compression as A, (A,)-,,, where 


4[0.357,(1 + g) + 1] 


(Ky)erit 
1.4gr,° + 1.307,(1 + g) + 1 


is the buckling load of the panel. We shall plot the 
stresses and deflection against this reduced edge com- 
pression A, (A,),;;: instead of the edge compression 
K,. In this way, many of the ten curves mentioned 
above will converge approximately to a single curve, 
which gives a decided advantage over the other way 
of plotting. 

The deflection ~z’,,,,, the membrane stress in the con- 
cave facing (7"y,) mar, and the shear stress in the core 
mar, are plotted nondimensionally against the re- 
duced edge compression in Figs. 3, 4, and 5. In each 
of these figures there are two sets of data, one repre- 
senting the effect of initial warping and the other the 
effect of transverse pressure. For each set of data there 
are four different subsets of points representing the 
parameters r, = 0, 0.1, 0.2, 0.3, for g = 1. The points 
that correspond to g = 0.8 and 1.25 are not shown be- 
cause in our plots these data come very close to the 
corresponding points for g = 1. All these curves show 
that, owing to the initial warping, the deflection and 
stresses start from zero at zero edge compression, as it 
should be, and increase rapidly with the increase of 
edge compression, and finally go to infinity when the 
edge compression approaches the theoretical buckling 
load. This situation would not actually occur, of 
course, because of the appearance of inelastic effects 
and because our small deflection assumption would no 
longer be valid. The deflection and stresses are not 
zero at zero edge compression because of the bending 
caused by the uniform transverse pressure, and they in- 
crease at a slower rate with the increase of edge com- 
pression than in the previous case. 

It is seen in Fig. 3 that, due to the effect of the initial 
warping, the plot of deflection vs. reduced edge com- 
pression is represented exactly by a single curve 


Ky, (Ky)erit 


A) = 
i= Ky, (Ky)crit 


for all the different panel parameters. Figures 4 and 
5 show, however, that the stress vs. reduced edge com- 
pression data does not fall on a single curve as does 
the deflection curve in Fig. 3. Rather, at a given 
value of the reduced edge compression, the larger the 
value of r, the smaller are the stresses. This is to be 
expected, since at a given value of K,/(K,),.,i: the de- 
flection of panels with different r, are the same, while 
panels with larger values of 7, experience greater shear- 
ing deformation and therefore carry smaller stresses. 
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As to the effect of the uniform transverse pressure, 
the situation is different. As shown in Fig. 3, the plot 
of deflection vs. reduced edge compression does not 
become a single curve. It is interesting to see, how- 
ever, that in Fig. 4, the plot of membrane stress vs. 
reduced edge compression for different values of the 
parameter r, almost exactly merges to a single curve, 
and in Fig. 5, the plot of shearing stress vs. reduced 
edge compression falls approximately on a single curve. 

The surprising result that the stresses due to the 
uniform transverse pressure could be represented by a 
single curve in the plot of stress vs. reduced edge com- 
pression is very significant indeed. It makes a great 
simplification in the stress calculations for sandwich 
panels. Since sandwich panels with different values 
of r, correspond to a single curve in the plot of stress 
vs. K,/(Ky)-ri, and in particular to the curve repre- 
sented by r, = 0, which corresponds to the usual thin 
plate, the numerous results obtained in the theory of 
thin plates can readily be applied to sandwich panels 
once we know the buckling load of a sandwich panel. 
The buckling loads of sandwich panels have been a 
subject of extensive investigation and are readily avail- 
able in the form of design charts.’ For the stresses 
caused by the initial deflection of the panel, Figs. 4 and 
5 show that the plots of stress vs. reduced edge com- 
pression for different values of the parameters r, do 
not fall on a single curve. However, these graphs do 
show that the stresses obtained for r, = 0 are upper 
bounds of stresses for other values of r, and therefore 
the results from the thin plate theory can again serve 
as a design criterion in the preliminary design of sand- 
wich panels. 
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Fic. 3. Maximum deflection coefficients A; and A» versus re- 
duced edge compression. 
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Fic. 4. Maximum facing membrane stress coefficients B,; and 
Bz versus reduced edge compression. 


The graphs also show that the ratio of the maximum 
shearing stress in the core to the maximum membrane 
stress in the facing is of the order of magnitude of 
mt”/L,. The linearly varying curvature stresses in 
the facings were computed, but are not shown here. 
It suffices to say that the results show that the maxi- 
mum values of these stresses is of the order of magni- 
tude of ¢”/7 times the maximum membrane stress in the 


facing. 
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Fic. 5. Maximum core shearing stress coefficients C; and C: 
versus reduced edge compression. 
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Discussion and Conclusion 


(a) The paper presents a practical method of solving 
the elastic problem of an initially warped sandwich 
panel under combined loadings. It can be applied to 
the simply supported sandwich panel under tempera- 
ture gradient, as long as the small deflection theory 
holds. The theory is so general that the sandwich 
facings may also be of different thicknesses and different 
materials. The sandwich core can be orthotropic and 
of very low shear rigidity. The loadings considered are 
biaxial compression and arbitrary transverse pressure. 
The method can be extended without difficulty to panels 
under eccentric edge compression and edge moments. 

(b) In this method the stresses and dedections of a 
sandwich panel are represented by doubly infinite 
series, or sums of a finite number of terms, depending 
on whether the transverse loading and initial warping 
of the panel are expanded as doubly infinite series or 
as sums of a finite number of terms. Therefore con- 
siderable reduction in computational work would be 
obtained if the initial warping and transverse loading 
could be approximated by just a finite number of terms 

(c) As an example, a square sandwich panel is 
chosen for numerical computation. The normal 
stresses in the facings are comparatively easy to calcu- 
late from the doubly summed infinite series. The 
convergence is good. The series representation for the 
facing shear stress converges very slowly. Fortunately 
such stresses are not important in design and need not 
be calculated. The doubly summed infinite series for 
these shearing stresses in the core converge more 
slowly and are more difficult to calculate. This com- 
putational difficulty is greatly reduced by our new 
method of presenting the data. It is shown that if 
the deflection and stresses are plotted against a reduced 
edge compression, defined as the actual edge com- 
pression divided by the critical edge compression of 
that panel, the data for panels with different physical 
parameters fall approximately on a single curve. 
Hence, a few computed curves can represent a large 
family of panels with different panel parameters r, 
and g. 


(d) With this new method of presenting data, the 
present calculations can be generalized to a large 
range of different physical parameters. The solution 
of sandwich panel with initial warping and _ trans- 
verse pressure can be obtained with ease. It is ex- 
pected that this simple presentation is especially use- 
ful to structural engineers for practical design purposes. 
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Derivation of an Analytical Expression 
for Supersonic Spillage of Single-Cone Inlets 


Craig W. Hartsell* 
Space Technology Laboratories, Inc., Los Angeles, Calif. 
May 9, 1960 


SYMBOLS 

A = area 

m = mass-flow rate 

y = radius measured from centerline of cone 

x = dimension measured along centerline of cone with x = 0 at the cone 
tip 

= tan 

6 = half angle measured between the centerline of the cone and a line 
emanating from the cone tip 

Subscripts 

c = cowl 

L = cowl lip 

Ss = cone surface 

o = free stream 

w = wave front 


SUPERSONIC capture-area ratio of single-cone inlets 
operating critically or supercritically has been found! for 
a limited range of Mach numbers and cone angles (Mo to 5, 
6, from 15° to 30°). It would often be useful to have an equa- 
tion which would provide the capture-area variation over a wide 
range of Mach numbers and cone angles. 

It has been shown? that the streamlines behind a conical shock 
may be approximated by portions of hyperbolas for a wide variety 
of conditions. With this assumption, the desired expression may 
readily be derived by making reference to Fig. 1. 

The streamline equation is 


* This work was performed while the author was connected with The 
Marquardt Corporation. 
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r? = 7,2x2 + constant (1) 
From Fig. 1, point (b) is the cowl lip and by using the cowl-lip 
coordinates the constant may be eliminated 
The capture stream-tube is bounded by the streamline extending 
from (b) to (a). At point (a), the coordinates define the cap- 
tured stream-tube, thus 


rot — ry? = — (3) 
Ao/A. = Aw/AL = rw?/rz? (4) 
tan 9 = rw/Xw (5) 
tan 0, = (6) 


Combining Egg. (3), (4), (5), and (6) 


Ao m cot? @, — cot? 0, 
= = (7) 
A, Wikis cot? 0, — cot? 6, 

A comparison between results from the use of Eq. (7) and 
Sibulkin’s! results are presented in Figs. 2 and 3 for two cone 
angles. It may be seen that the present approximation yields 
an excellent correlation. As the flight Mach number increases, 
the correlation improves until at moderate supersonic speeds 
identical results are obtained. The usefulness of this approxi- 
mate solution is considered to be restricted as follows: 


Fic. 1. Geometry of single-cone inlet. 
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6, > 10°, My > 1.5 
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A Comparison of Blunt-Body Flow-Field 
Results} 


Herbert W. Ridyard 
Aeronautical Research Engineer, Missile and Space Vehicle 
Department, General Electric Co., Philadelphia, Pa. 


April 11, 1960 


INTRODUCTION 

N A RATHER COMPREHENSIVE STUDY of the supersonic blunt- 

body problem, Van Dyke! presented a numerical method 
which determines the flow field for a family of plane and axisym- 
metric detached shock shapes. Van Dyke's results! were shown 
to agree quite well with experiment and with the numerical re- 
sults of Garabedian and Lieberstein.?, Recently, Vaglio-Laurin 
and Ferri*® introduced another numerical solution for the blunt- 
body problem, and their results in the subsonic region were 
stated to compare satisfactorily with the data of reference 1. 

In order to continue the blunt-body solution further down- 
stream in the supersonic region, the results of a forward-integra- 
tion method in the neighborhood of the sonic line may be used as 
input for the method of characteristics. It is, therefore, of some 
interest to compare these methods in this region of the flow. This 
note presents some of the results from such a comparison between 
the Van Dyke forward-integration method and the method of 
characteristics. 


APPLICATION OF THE METHODS 
Both the Van Dyke and the characteristic calculations were 
performed at a Mach number of 20 assuming an ideal gas (ratio 
of specific heats = 1.4) for a given parabolic shock shape defined 


+ This work was performed under the auspices of the USAF Ballistic 
Missile Division, Contract No. AF 04(647)-269. 
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by the relation, r? = 2Rsx where r and x are the radial and axial 
components, respectively, of a cylindrical coordinate system and 
Rg is the radius of curvature of the shock at the axis. The extent 
of the flow-field calculations is illustrated in Fig. 1. The entire 
mixed-flow region, bounded by the solid lines A BCD, was deter- 
mined by the Van Dyke method. The characteristic triangle 
A’B'P’ in the supersonic region (Fig. 1), bounded by the shock 
and the dashed characteristic lines, was obtained essentially by 
the procedure outlined in reference 4. 

Following the method of reference 1, the shock-oriented mesh 
for the Van Dyke calculation had fifty-five steps along the shock 
in the £ direction (see Fig. 1) with AE = 0.05. The computation 
required six steps in the 7 direction normal to the shock (An = 
0.025) to determine the body. An 11-point differentiation for- 
mula with central differences was used to determine the £-deriva- 
tives of the flow quantities at each value of n. In order to avoid 
other than central differences, five points were dropped at the 
downstream end of the mesh at each step in the yn direction (as 
suggested in reference 1). Consequently, only 20 points were 
obtained at the body, curve CD, Fig. 1; the corresponding 20 
points at the shock define the curve AB. 

The characteristic net A’B’P’ (Fig. 1) was determined by use 
of 60 points along the segment of the shock A’B’ in even incre- 
ments of the nondimensional radial coordinate, r/Rs, with 
A(r/Rs) = 0.005. One iteration of the solution was performed 
at every point in the net. 


COMPARISON OF RESULTS 
Since the Van Dyke and the characteristic calculations were 
performed for identical shock shapes and free-stream conditions, 
the flow quantities immediately behind the shock were precisely 
the same for either flow field. It is, therefore, of interest to com- 
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nondimensional entropy as a function of y/é for the method of 
characteristics and the Van Dyke method. 


TABLE 1 


Angle 
Mach No. (deg.) 


Method x/Rs r/Rs 


Characteristics 


at body point P’ 0.4040 0.6466 1.079 45.33 
Van Dyke at 

body point P 0.4080 0.6447 1.076 44.51 
Per cent difference 1.0 0.3 0.3 1.8 


pare the two solutions as they proceed away from the shock. 
Such a comparison has been made by use of interpolated flow- 
field results along a line OP drawn from the shock at point O to 
the ‘‘Van Dyke’”’ body at point P in the y direction indicated 
in Fig. 1. This line OP, which includes the body point P’ of the 
characteristic calculation, was constructed normal to the flow 
direction at P’. The results of the comparison are indicated in 
Fig. 2 for three of the flow quantities—Mach number, M, flow 
angle, 0, (measured with respect to the axis of symmetry) and the 
nondimensional entropy S/R as a function of the nondimensional 
distance from the shock, y/6._ Here R is the gas constant and 6 
is the total distance from the shock to the ‘‘Van Dyke’ body 
along the line OP in Fig. 1. 

As shown in Fig. 2, there are small differences between the re- 
sults which increase in magnitude as the distance from the shock, 
y/6, increases. The values of the entropy S/R for the two 
methods at the respective body points are, of course, identical 
since they must coincide with the stagnation-streamline boundary 
condition. 

A detailed comparison of the results at the body points P and 
P’ (where the maximum differences occur) for the two methods is 
presented in Table 1. The per cent differences are seen to be small 
except perhaps for the flow angle 6, for which the per cent dif- 
ference is about 1.8 per cent (A@ = 0.82°). 


CONCLUSION 

In general, the present comparison indicates that the Van Dyke 
forward-integration method is in good agreement with the con- 
ventional method of characteristics in the supersonic-flow region 
in the neighborhood of the sonic line. This is a rather interesting 
result since the Van Dyke calculation was performed with only 
six integration steps as compared to about sixty for the charac- 
teristic calculation. It must be noted, however, that no attempt 
was made to investigate the effect of the integration step size in 
this exploratory investigation. 
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A Comparison of Transient and Quasi-Steady 
Performance of Melting-Type Re-Entry Shields 


Ernst Wilhelm Adams 
Aeroballistics Laboratory, ABMA, Redstone Arsenal, Ala. 
June 10, 1960 


eo INVESTIGATIONS of the ablation-type heat- 
protection problem usually are restricted to a small vicinity 
of the stagnation point, because there the dependency of the 
problem on the coordinate parallel to the surface can be elimi- 
nated by employing symmetry considerations. The published 
investigations of the said problem present quasi-steady solutions 
only, which can by obtained—e.g., by use of the references 2 
and 3. In reference 1, an exact transient-solution procedure for 
the problem at hand has been published recently. The present 
note gives a comparision of evaluations of reference 1 and refer- 
ence 3 for both a ballistic satellite and an IRBM, which are 
protected by homogeneous, opaque, and nonevaporating melting- 
type glass shields. The material properties of these shields are 
given in Table 2. The trajectory data of the vehicles are pre- 
sented in Table 1 and in Fig. 1, which shows the flight speeds V 
(dash-dotted lines) as functions of the altitude H and the time ¢. 

The Figs. 1-3 present some performance data, which are ex- 
plained in the legends, as functions of the altitude H. Fairly 
good coincidence of the quasi-steady and the transient solutions 
takes place in the vicinity of the maximum temperature only. 
This is explained as follows: At every point in the shield, the 
heat conduction and the heat convection (due to the motion of 
the glass particles relative to the surface) cancel each other under 
steady conditions. By definition, heat transfer from the shield 
into the air cannot take place in the steady state. Hence, the 


t (SEC) "SATELLITE" 


t_[sec] "iRBM 


Fic. 1. Flight speed V (dash-dotted lines) and surface tempera- 
ture 7;. 


Dashed lines: Quasi-steady solution for 7;. 
Solid lines: Transient solution for 7;. 
Heavy lines: Satellite results. 

Light lines: IRBM results. 
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Fic. 2. Melting speed -v7.. The meaning of the solid 
and dashed lines is explained in the legend to Fig. 1. 
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Fic. 3. Distance a from the surface to the station where the 
temperature is 811°K. = 1000°F. The meaning of the solid 
and dashed lines is explained in the legend to Fig. 1. 


quasi-steady solutions are poor approximations to the transient 
solutions during the initial heating period of the shield before 
melting starts, and during the final period after the transient 
heat flux across the surface has changed its direction. 

Melting of the shields takes place at and close to the surface 
and, hence, is governed by the surface temperature. The quasi- 
steady solutions considerably overestimate the total ablation 


s= — Sfv.(t) dt 


which is included in Table 1. 


TABLE 1 
Dimension Satellite IRBM 

Nose radius m. 1 1 
Ballistic factor m.3/kg. sec.” 30 1073 3.5 % 10-3 
Re-entry angle 92 .35° 124.9° 
Total ablation 

(transient ) mm. 15.5 3:3 
Total ablation 

(steady ) mm, 23.9 4.4 


keal./m. °K. sec. 
keal./kg. °K. 


k = 9.713 X 1074 


Cp = 0.29 


Thermal conductivity 
Specific heat 


Specific weight y = 2,100 kg./m.3 
Emissivity constant e =08 Re 
Vapor pressure po = 0 kg./m.? 


kg. sec./m.? 
= 0.017 X%exp {(4,230/T°K. )?-%64} 


Viscosity function 
of Pyrex glass 
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After the end of the initial heating period of the shield, the 
thermal state at its surface is determined by the strongly non- 
linear temperature dependency of the heat-transfer rates across 
the surface. The deviations of the quasi-steady and the transient 
solutions increase with the distance from the surface, as is seen 
by comparision of Figs. 1 and 3. This shows that the quasi- 
steady solution considerably overestimates the thermal pene- 
tration of the shield. 

The comparision with a pertinent transient solution indicates 
that a quasi-steady solution of the problem at hand yields un- 
reliable results and, in particular, fails regarding the necessary 
thickness of the heat shield—i.e., the sum of its total ablation 
and its thermal penetration at impact time. It is to be expected 
that similar conclusions hold for different types of re-entry 
vehicles. 
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Bending Stresses Due to Temperature 
in Hollow Circular Plates, Part | 
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T# PURPOSE OF THIS SERIES OF NOTES (Parts I-III) is to pre- 
sent formulas and curves for the deflection, moments, and 
shears in a circular plate with a concentric hole subject to a linear 
thermal gradient through the thickness. The various combina- 
tions of conventional boundary conditions at the inner and outer 
radii will be considered. The boundary conditions to be con- 
sidered in Part I are the cases of clamped outer radius coupled 
with clamped, simple, or free inner boundary. 
A diametrical section of the plate is shown in Fig. 1. The 
equilibrium deflection equation for an elastic plate subject to a 
linear gradient through the thickness takes the form (reference 1) 


Viw = 0 (1) 
where V‘ for the axisymmetric case is given by 
V4 = [(d?/dr?) + (1/r)(d/dr)] [(d?/dr?) + (1/r)(d/dr)| 
The general solution of Eq. (1) is 
w= + Cr? + C3 In (r/a) + Cyr? In (r/a) (2) 


where the coefficients C; are determined from the appropriate 
boundary conditions. 

The radial and tangential moments per unit of length J/, and 
Meare determined from! 


dw rdw (1+ 


M, = —D : 
dr? r dr h Tp (3) 
1 dw d*w (1+ 
Me = —D\| - — — — —— 
_ yr dr aiid dr? h Tp (4) 


where 7p is the temperature difference between the upper and 
lower faces respectively, w is positive downward, M, and Mo are 
positive when the upper fibers are compressed, v is Poisson’s ratio, 
and a is the linear coefficient of expansion. The radial shearing 
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Fic. 1. Nondimensional deflection and bending moments 
for plate clamped at the outer radius and free at the inner 
radius. 
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moments for plate clamped at the outer radius and simply sup- 
ported at the inner radius. 
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force per unit of length is given by (tangential component is zero 


from symmetry ) 
djid dw 
aD — 5 
dr E dr )] (5) 


where Q is positive when the resultant force on an elementary ring 
is downward from the inside and upward from the outside. 
Solutions are now given for the above edge conditions. 


Case A: Plate Clamped at Outer and Inner Radii 


w= Q=0, M,/EaTph? = 
Mo/EaT ph? = 1/{12(1 — v)] (6) 
Case B: Plate Clamped at Outer Radius and Free at Inner Radius 
The nondimensional deflection, moments, and shear are 

wh —(1 + »){1/2[1 — (r*/a*)] + In (r/a)} (7 
aTpb? {1 — v) + (1 + v)(b?/a?)] 
M, = (b?/r?)] (8) 

EaT ph? 12[(1 — +) + (1 + »)(b?/a?)] 
Me + (b?/r*)] (9) 

EaT ph? 12((1 — v) + (1 + v)(b?/a?)] ; 
Q=0 (10) 


Fig. 1 presents these quantities (multiplied by suitable powers 
of 10) versus r/a with b/a = 0.2, 0.4, 0.6, 0.8. 


Case C: Plate Clamped at Outer Radius and Simply Supported 
at Inner Radius 

The nondimensional deflection, moments, and shear are given 
by 


wh 1 a? 41 r 4 
aT pb? b? a’ b 
b r 
(11) 
b? a a 
M, 1 2 1 b? 491 b + 
=— — ——— 2 in | in - 
EaT ph? 12/11 — v)A a? a b 
b? b a? b? a? 


Mo af 
EaTph? 121 — »)A a? 


1+v/L \r a r2 a? |f 
Qb 1 b “| b 
= ———_} 1 21 ——j|- 14 
EaT ph? 3(1 — v?2)A a a? |r (14) 


where 


it. 


The results of Eqs. (11)-(15) are shown in Fig. 2. It should 
be remarked that, when the plate is subjected to both normal 
loads and temperature, then, for small deflections, the principle 
of superposition is applicable.? 


REFERENCES 


1 Timoshenko, S., Theory of Plates and Shells, 1st Ed., McGraw-Hill Book 
Company, Inc.,1940. 

2 Zaid, M., and Forray, M., Deformation and Moments in Elastically Re- 
strained Circular Plates Under Arbitrary Load or Linear Thermal Gradient, 
J. Mech. Eng, ASME Trans., No. 59-A-39 (Preprint), to be published in 
1960. 


FORUM 793 


Variable-Thrust Rocket Engines and Their 
Modes of Operation 


Mario William Cardullo 


U.S. Naval Air Rocket Test Station 
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SYMBOLS 


A = area (in.?) 

Ca = orifice coefficient 

c* = characteristic velocity (ft. per sec.) 

Cy = thrust coefficient 

F = thrust (Ib.) 

Is) = specific impulse (sec.) 

k = ratio of specific heats 

P = pressure (psi) 

S = bulk specific gravity 

V = injection velocity 

w = total flow rate (Ib./sec.) bulk density (Ib./in.*) 
Subscripts 

c = chamber 

e = exit 

t = injection 

t = throat 


INTRODUCTION 
a the past few years many types of controls have been devel- 
oped to give liquid-propellant rocket engines a wide range of 
operation. The purpose of this note is to analyze the various 
methods of achieving variable operation of liquid rocket engines. 
At present, the following methods are used to vary thrust: 
(3) Throat throttling 
(4) Combination of throat and 
injection throttling 


(1) Stepped chambers 
(2) Propellant throttling 


All of these methods have been tried and proved successful, to 
some degree. The performances of these devices vary con- 
siderably. 
ANALYSIS 
The combustion-chamber pressure (?,) for a rocket engine! is 
given by 
a 
P, = we*/Arg (1) 


where w can be found from 


w= CipiVA; (2) 

and V is given by 
V = (24 gAP/p;)'* (3) 

AP is defined as 

AP = P; — P. 
Letting 

a = 24C,?gp;, B = ac*/g 

and 56 = aP; 


Substituting Eqs. (1) and (2) into (3) yields the quadratic 
equation 
w? + — 6A? = 0 (4) 
Solving Eq. (4) 
w = (B/2)A;{ [(Ai/At)? + Nopro.]'/2 — (Ai/Av)} (5) 
where: Noro. = 304 P;/S c** 

The propellant number (Npro.) is based upon a nominal value 
of C, of 0.7. Table 1 shows values of P;/Nopro. for various com- 
mon propellant combinations. 

Substituting for w into Eq. (1) 


P./P; = (2/Npro. \(A;/A,){ + (A,/A:)} 
(6) 


Eq. (6) has been plotted in Fig. 1 for various propellant num- 
bers. It can easily be seen that unless the ratio (A;/A;) remains 
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Values of P;/Npro. for Various Propellants* 


Oxidizer Fuel Noro. (psi) 
IRFNA UDMH 1.24 X 10° 
H202 (90%) RP—1 
F, 1.04 105 
N2O, NoH, 1.11 X 10° 
Oz C:H;OH (92.5%) 1.03 X 105 


* Based upon shifting equilibrium values from references 
(2) and (3). 


TABLE 2 
Results of Throttling 
Engine Parameters Specific Impulse 


| Pressure Vv 


Method | 

| Constants Variables acuum |Atmospheric 
Pi, Ae, |Various 
| At, engines 


Stepped engines Constant) Constant |Constant 
| 


Propellant |Ai,Ac, Ae |Pi Decrease Constant Decrease 
|Pi, Ae, At |Ai Decrease | Constant | Decrease 
Throat throttling |P;,Ai,Ae (|At Increase | Increase |Increase 


Throat and 


A-/Ai, |Ae/Ai Constant) Increase |Increase 
injection As, Pi | then 


throttling | | decrease* 


* Depends upon initial design conditions. 


constant the chamber pressure will change as the engine is being 
throttled. 
The thrust of a rocket engine is given by 


F = C,w (c*/g) (7) 
where C, is defined by 
Cy = [(2k*/k — 1)(2/k + + W/E 
[1 — Vk + (P, — 


The ratio P,/P, is purely a function of A;/A, and thermody- 
namic properties of gases. Since A;/A, is an implicit function of 
the pressure ratio, we can only express the thrust in the functional 
form 


F = f(A1/Ai, A,/A%) 


The specific impulse can also be placed in a similar form. 
Fig. 2 shows the effect of the throttling parameters of the 
specific impulse for vacuum conditions. 


DISCUSSION 


The various devices originally mentioned can now be analyzed 
on the basis of the equations developed. Table 2 shows, quali- 
tatively, the result of throttling on chamber pressure and spe- 
cific impulse for the various types of variable-thrust devices. 


CONCLUSIONS 


The only device that will maintain constant chamber pressure 
under all throttling conditions is that engine where the ratio 
A,/A; is designed to remain constant. The only limitation of 
this system is that imposed by expansion conditions. This may 
possibly be remedied by incorporating a ‘‘plug’’ nozzle in the 
design. Thus the ratio A-/A; will also be constant, where A, is 
exit flow area for a ‘‘plug’’ nozzle. This engine will have con- 
stant specific impulse under all operating conditions and will be 
capable of thrust variations of any degree. 
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T# AIRLOADS ON WINGS due to strong blast waves have fre- 

quently been estimated by linearized methods of the type 
employed in gust analyses. The need for experimental support 
of this procedure has recently been emphasized by Ehlers and 
Shoemaker,' and Goodman and Sargent.?, A program of meas- 
urement of the transient loads on airfoils and wings exposed to 
simulated strong blast waves in a shock tube has been underway 
at M.I.T. for several years,* and experimental results have been 
published in reference 3. In reference 4 these measurements 
have been compared with the traveling-gust theory of Miles,5 
Hobbs, and Drischler and Diederich.’ The purpose of this note 
is to review several of the principal results of the latter com- 
parison. 

The application of the traveling-gust theory to blast waves 
was proposed independently by Hobbs,® and Drischler and 
Diederich’ by equating the velocity of the traveling-gust front to 
the velocity of the shock fronting the blast wave. The gust- 
induced angle of attack is equated to the change in angle of 
attack experienced by the wing as it penetrates the shock front. 

The traveling-gust model appears to be a reasonable representa- 
tion of a weak blast wave. But since it is based on the acoustic 


* This investigation has been sponsored by the Aircraft Laboratory of the 
Wright Air Development Division. 
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equations of motion, there is serious question whether the theory 
can be adapted to strong blast waves. In the reported tests, the 
shock waves were strong: the Mach number of the shock- 
induced flow M2 was (.4-1.0, and the shock-induced angle of 
attack @ was as high as 30°. Therefore, both the departures of 
the loading in the vicinity of the strong shocks, and the effect of a 
large a were investigated. 

The chordwise distribution of the load on an NACA 65-010 
profile while the shock front is diffracted by the airfoil is plotted 
in Fig. 1. The loading coefficient Ac, is equal to Ap/(1/2)p.u2?, 
where Ap is the pressure difference across the airfoil, pz and uw» are, 
respectively, the density and material velocity behind the shock, 
x is the chordwise distance from the leading edge, a» is the speed 
of sound behind the shock, and ¢ is the time after shock impinge- 
ment. For the traveling-gust theory Ac, was computed using 
two different coordinate systems: (1) the coordinates were 
fixed on the fluid ahead of the shock, and (2) the coordinates were 
fixed on the free stream behind the shock, respectively designated 


Corg—1 and Cpig-2 The lift is plotted in Fig. 2 for two profiles, 
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where c; = (section lift)/(1/2)p2u2*c, and ¢ is the chord. The 
functions Chey and correspond to the traveling-gust theory 
with respective coordinate systems (1) and (2). In both figures 
the functions have been reduced by a. 

The first conclusion is that the shock-induced coordinate sys- 
tem (2) provides considerably better agreement with both the 
load distribution and the lift than the initial flow coordinate 
system (1). Secondly, both Ac,/a and c;/a@ are essentially inde- 
pendent of a while the shock is bifurcated by the airfoil (0 < 
ust/¢ < 0.3), and in good agreement with the theory in the shock- 
induced coordinates, with the exception of the load distribution 
near the shock front. (The loading functions were computed for 
regions not affected by the trailing edge.) The influence of com- 
pressibility, or specifically the speed of sound, is demonstrated by 
the curve designated 


M = 0)/V1 Me? 


where Ctygl M = 0)is Clips for M, = 0. It is apparent the com- 
pressibility cannot be handled during this early time by the 
Prandtl-Glauert factor. Finally, the effect of the gust speed is 
demonstrated by the two curves designated Ct, and ¢;,, which 
The 
data indicate that the effect of the gust speed on the loading 
became small by w#3f/c = 1. Thereafter the influence of viscosity 
is important at large angles of attack,’ which is the subject of 
current investigation. 


represent ¢7,)_, with a zero and infinite gust-front velocity. 
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Comments on ‘‘Mach-Line Determination 
for Air in Dissociation Equilibrium’”’ 


Eldon L. Knuth* 
Research Consultant, The ASTRO Division, 
The Marquardt Corporation, Van Nuys, Calif. 


April 18, 1960 


Keo PURPOSE OF THIS NOTE is to clarify and to supplement 
the results presented in reference 1, in which note the condi- 


+ 


tions necessary for constant-density flows through shock waves 
are examined. In the present note it is shown that these condi- 
tions do not lead, in general, to Mach waves. 

Since Mach waves are weak disturbances, it is necessary to 
examine the case in which small variations exist across the wave. 
Substitute, therefore, 


p2/p. = 1+ dp/p 
K2/K, 1+ dK/K 


into the expression for the density ratio across a shock! 


p2 2 


to obtain, neglecting terms of second and higher orders in differentials 


2 


9 


2 


2 (1 — K) — (1+ JK 


* Also, Associate Professor of Engineering, University of California, Los Angeles, Calif. 


The author wishes to express his appreciation for the stimulating technical discussions with W. J. Hayes of the ASTRO Division, The Marquardt Corp. 
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with 

p = density 

y = isentropic exponent = (01n p/d1n p)s 

\J_ = Mach number based on velocity component normal to 
wave 

K = (p/p)/(h + e) 

p = pressure 

h = enthalpy 

é = internal energy 


and with the subscripts 1 and 2 referring, respectively, to condi- 
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disturbances may be treated as isentropic flows, one may write 

dK dp dp dh+e) 

Kp p hte 
_dp dp _ p/p **) 
p p h+t+e\p p 

= — 1) — K(y+ 

Substituting into the linearized equation, one obtains, using the 
upper signs, 


(1 — de 


tions before and after the wave. Since flows through weak : — =@ (1) 
(i — — (1 + 
and, using the lower signs, 
p 
—1) — K(y + 1)] | K + yM2) + 
yM,? + [(y ) (y + 


These two equations correspond to the two roots of the quadratic 
density-ratio equation for the case in which small variations 
exist across the wave. 

Eq. (1) corresponds to Mach waves. 
tion is satisfied if 


For dp # 0, this equa- 


M;?? = 1 (3) 


in agreement with the definition of a Mach wave. 

Eq. (2) describes the density change across weak shock waves. 
For perfect gases, neglecting terms of second and higher orders 
in (7,? — 1), this equation reduces to 

dp/p = 2(M,? — 1)/(y + 1) (4) 
in agreement with results found in gasdynamics texts. 

For real gases, if no density change occurs across the shock 
wave, one obtains, from Eq. (2), the result given in Eq. (7) of 
reference 1—i.e., 


M,? = h/(ye) (5) 
which, for perfect gases, reduces to 
M;2=1 (6) 


Eq. (5) indicates that, for real gases, constant-density flows 
through shock waves are possible for 4, # 1. Consequently, 
these shock waves with constant-density flows are not, in gen- 
eral, Mach waves. 
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Remarks on the Eddy Viscosity in Compressible 
Mixing Flows 


Lu Ting* and Paul A. Libby** 


Polytechnic Institute of Brooklyn, and 
General Applied Science Laboratories, Inc. 


May 9, 1960 


I CONNECTION WITH A STUDY of the wakes behind bodies in 
hypersonic flow carried out for the Missile and Space Vehicle 
Division of the General Electric Company, it was desired to 
estimate the eddy viscosity in axisymmetric, compressible 
wakes. Because of the lack of applicable experimental data, 
it was found necessary to make such an estimate by rationally 
extending the few available data for incompressible flows to the 
compressible case. This suggested the application and ex- 

* Research Associate Professor and Scientific Investigator. 

** Professor and Senior Scientific Investigator. 


tension of the transformations applied to turbulent boundary 
layers in reference 1. Making the assumption that the turbu- 
lent shear stresses over an infinitesimal mass are invariant with 
transformation, Mager! showed that the partial differential equa- 
tions for the compressible turbulent boundary layer can be 
transformed to incompressible form. The validity of this as- 
sumption and of the transformations was established for several 
boundary-layer flows by comparison with experiment. 

Consider the application of these concepts to two-dimensional 
mixing flows, e.g., wakes or jets; in this note the symbols of refer- 
ence 1 are used for convenience. It is desired to find a relation 
between the compressible and incompressible eddy viscosities 
denoted by e«, and e,*, respectively; it follows from Mager’s 
analysis that: 

Tw = pe(Ou/dy) = pu'v’ = F*pou'y’* = F*e,*(Ou*/OV) (1) 
With u = Fu* and 0/Oy = F(p/po)d/dy there results 

(po/p)?e:* (2) 


where po is a reference density for the flow, e.g., the density in 
one stream external to the mixing region. It should be noted 
that the scale factor G introduced by Mager is unnecessary for 
wake flows if, as is customary, the laminar shear is neglected. 

Now, it is generally assumed in theoretical analyses and con- 
sidered to be confirmed by experiments (cf. reference 2) that the 
incompressible eddy viscosity ¢,* is constant across the mixing 
zone, i.e., independent of the transverse coordinate Y. The 
appearance of (po/p)? in Eq. (2) implies that the compressible 
eddy viscosity «, on the other hand is not independent of the 
transverse coordinate y. Thus, application of Mager’s analysis 
to two-dimensional mixing problems leads to a prediction con- 
cerning the compressible eddy viscosity in disagreement with 
that usually assumed for mathematical convenience in analyses 
of compressible mixing problems (cf. references 3 and 4). 

Consider next the application of Mager’s transformations to 
axisymmetric wakes. As might be expected on physical grounds 
the assumption concerning the invariance of the turbulent shear 
stress over an infinitesimal mass must in this case be altered; 
the corresponding assumption is that the moment about the 
axis of the turbulent shear stress over an infinitesimal mass is 
preserved; i.e., 


(rpu'v’)p(r dr dx) = (r*pou'v’*)( por*dr*dx*) (3) 


With this assumption the partial differential equations for an 
axisymmetric wake of a compressible fluid are transformed to 
incompressible form, and the eddy viscosities are related by the 
equation 


€1*( po/p)*(r*/r)? (4) 


r 
where r*2 - ff 2( p/po)rdr 
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The validity of both Eqs. (2) and (4) must be established by 
experiment and must, therefore, be considered provisional. 
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The Flow Field in the Diffuser 
of a Radial Compressor 


Inge L. Ryhming* 
Fluid Dynamics Research Group, Convair, San Diego, Calif. 
May 9, 1960 


beg NOTE discusses the two-dimensional diffuser flow field 
in a radial compressor outside the impeller wheel. It is 
assumed that the diffuser has guide vanes arranged in a circular 
row at a radius 72. The impeller wheel has the radius 7, (see 
Fig. 1). The flow in the diffuser starts at the circle with the 
radius 7;. The velocity components, u; and 1, in the r and 6 
directions of the velocity vector q; on this circle are prescribed 
together with the thermal state of the gas The flow so prescribed 
on the radius 7, will, if no disturbances are present (i.e., no 
boundary conditions in the flow other than zero velocity at in- 
finity are to be fulfilled), develop in a spiral flow. 

The flow equation for two-dimensional compressible flow 
written in polar coordinates, 7 and 8, is 


u,(a? — + = — 
(1/r)uvity — uvv, + a%(u/r) = 0 (1) 
where a is the speed of sound. A solution independent of @ and 


satisfying the given initial conditions on the radius 7 is desired. 
Hence with 0/00 = 0 for irrotational flow we have 


ur(a? — u?) — uvv, + a2(u/r) = 0 (2) 
+ v = (or), = (3) 


From Eq. (3) it follows that 
vr = const. = yr, = A (4) 
The speed of sound a is given by the energy equation 
+ 9? + [2/(y — = [(y + 1)/(y — 1)Ja*® (5) 


where y¥ is the ratio of the specific heats and a* the critical speed. 
Eliminating a and v between the Eqs. (2), (4), and (5), an 


* Formerly, Staff Member, Flight Sciences Laboratory, Boeing Scientifi- 
Research Laboratories. 


Fic. 1. 


Arrangement of the diffuser. 


¥+7=1.52 

or M=1.30 
v+7=1.34 
or M=1.20 


Fic. 2. Integral curves of Eq. (7). 


Fic. 3. Flow field in the diffuser, vo(7,) = 0.2 and = 1.32. 


ordinary differential equation is obtained for u = u(r). By 
introducing new dependent and independent variables in this 
equation 

vy = (u/a*)?, n = (A/a*r)? (6) 
Taylor’s! equation results: 
& — ie — Be 


dy ny +1 —(y—1)n 


(7) 


Typical integral curves in the v, 7 plane of Eq. (7) applicable to 
the present problem are shown in Fig. 2. Along the straight 
line defined by » + 7 = 1 the Mach number of the flow is equal 
toone. To the right of this line in the v, 7 plane the Mach num- 
ber of the flow is greater than one. 

In the diffuser, however, we have assumed that guide vanes 
are present. We make the further assumption that the thick- 
ness of these guide vanes and their angle of attack with respect 
to the oncoming stream are small, so that only small disturbances 
are introduced in the spiral flow as described above. 

Suppose that the Mach number of the flow on the radius 1; 
is greater than one. Then a number of different possible flow 
types are possible. Since the flow Mach number decreases 
with increasing radius, the flow regime between the impeller 
wheel and the guide-vane cascade can be either fully supersonic, 
or, partly supersonic and partly subsonic. In the latter case, the 
disturbances from the guide-vane cascade will influence the entire 
flow field. This case will not be further considered here. 

In the case where the flow is fully supersonic, a series of waves 
starting from the leading edges of the guide vanes will develop. 
Depending on the Mach number and the flow direction at the 
radius 72, the spiral flow between the radii 7; and 72 can be either 
disturbed or undisturbed. In the latter case the waves are im- 
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mediately reflected inside the guide-vane cascade, whereas in the 
former case a wave from the leading edge of any guide vane can 
travel ahead of the leading edge of the next guide vane in the 
cascade. 

In the literature, a number of authors” * have considered this 
problem for a two-dimensional straight cascade immersed in a 
uniform semi-infinite parallel supersonic stream. In the present 
problem, however, the oncoming flow, i.e., the spiral flow, is 
nonuniform and of finite extent. Therefore, a flow with waves 
originating from the guide vanes between the radii 7; and r2 can 
only be possible if the initial data on the radius 7; are compatible 
with the boundary conditions on the guide vanes. That is, 
the initial data at the radius r; cannot be chosen arbitrarily since 
the flow at 7 is influenced by the boundary conditions at rz. 
Because of the fact that the oncoming flow is of finite extent, the 
waves eventually reach the radius 7, and are reflected. After their 
reflection at the radius 7; they may travel back to the radius r2 
and be reflected inside the guide-vane cascade. The existence 
of such a flow model is essentially dependent on the boundary 
conditions on the guide vanes. 

The reflection at 7 is physically complicated since the waves 
are reflected at moving boundaries, that is, on the rotating im- 
peller-wheel cascade. Unstationary effects of a type similar to 
those discussed in reference 4 will, therefore, occur in the flow 
field in this case. 

For a closer examination of the flow field, Eq. (1) has to be 
used. With the assumption of small disturbances, we can write 


+ iu, + 3, a=atda (8) 


where the quantities with the zero subscript satisfy the spiral- 
flow solution, Eqs. (7), (4), and (5). Introducing the Eqs. (8) 
in Eqs. (1) and (5), and retaining first-order terms in @ and 3 


gives 
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Solutions can now be found by the method of characteristics. 
Because the coefficients of the unknown terms in Eq. (9) are 
known functions, the characteristics are known in the flow field. 
A typical flow pattern is shown in Fig. 3. 

The flow field as discussed above can break down, however, 
if the Mach number at the radius 72 is too close to one, so that 
transonic effects occur, or if the shock waves go through several 
reflections between the impeller wheel and the guide-vane cas- 


cade. 
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An Approximate Investigation of the Effect 
of Boundary-Layer-Control Pumping 
on Powerplant Performance 


T. F. Kirkwood 
Aeronautical Engineer, The RAND Corporation, Santa Monica, Calif. 


May 17, 1960 


INTEREST in the application of boundary-layer 
control to aircraft has prompted the first-order analysis 
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of the effects of boundary-layer-bleed pumping on the perform- 
ance of turbojet and turboprop powerplants described here. It 
has become conventional! to define an ‘equivalent pumping- 
drag coefficient”’ as 

Co, = D./qaS (1) 


where D, is the ‘‘suction drag’’—i.e., the loss in thrust from the 
powerplant due to the diversion of power required to run the 
boundary-layer pump. 

For the purposes of this analysis the pumping system is as- 
sumed to consist of a pump driven either from the main turbine 
shaft, or, from an auxiliary turbine located downstream of the 
main turbine and having the same efficiency. Thrust is ob- 
tained from the pump exhaust, and the pump-pressure ratio is 
optimized to give the least net thrust loss. 

The suction drag is evaluated from the following considera- 
tions. The amount of power diverted to run the pump is equal 
to the work done on the pumped air per unit time divided by the 
efficiency of the pump: 

pump power = cp,ATW4/ns (2) 

If all of this power were completely lost, the suction drag would 
be given in terms of the engine propulsive efficiency », and flight 
speed as 

D, = (CpATWo/ )(np/ns) (3) 
However, there is the possibility of obtaining some thrust from 
the pump exhaust. If the pump exhaust velocity is U, the 
thrust gained in this way is W,/g(U — U,,) and the net suction 
drag becomes 
D, = [(cpATWo/U «)(np/ns)] [(Wa/g)(U = Ve )] (4) 
Introducing a pumped-airflow coefficient defined as 
Cw = We/pogsUa (5) 
and expressing the suction drag in coefficient form results in 
Co, = — [((2Cu/U.)(U — (6) 

The pump must take air from its static pressure at the surface 
of the wing, overcome a pressure drop in the slot and in the duct 
leading to the pump, raise the air to ambient pressure, and then 
raise it to a higher pressure sufficient to result in the desired 
exhaust velocity U. The temperature rise required for these 
purposes can be expressed approximately in terms of duct-pres- 
sure recovery 7, and the average wing-pressure coefficient 
Ap/qas 

U? 
Af a + 7, + 


2 


-1 y¥-1 A 
[= wat (2*)] (7) 
q 


The expression for the suction-drag coefficient given in Eq. (6) 


then becomes 


Np U 2 ( ) 
Mo? 


1-1 U 
(22: 2c. ( — :) (8) 
q 2 Ve 


The pump-discharge velocity U is optimized to give the mini- 
mum suction-drag coefficient by differentiating Eq. (8). The 
optimum discharge velocity is found to be given by 


Vont./U ns/Np (9) 


This relation implies that, for turboprop engines, where 7p 
and 7, are approximately equal, the optimum pump-discharge 
velocity is equal to the free-stream velocity. For turbojet 
engines at low flight speeds, the propulsive efficiency will be less 
than the pump efficiency, and Eq. (9) indicates that a discharge 
velocity greater than the flight velocity should be used. 

If the optimum discharge velocity is substituted into Eq. (8), 
the minimum suction-drag coefficient is found to be 
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Pressure loss 
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Fic. 1. Breakdown of suction-drag coefficient as a function of 
flight speed and engine type. 


The terms of Eq. (10) have the following significance. 
Momentum Increment: 
2 — 


This represents the change in thrust due to difference in propulsive 
efficiency of the main engine and the pump-discharge jet. On 
turboprop engines, where 7, ~ ns, this term is one. In this situa- 
tion, the pumped air is discharged at flight velocity and this 
term represents the loss associated with raising the total pressure 
of air from ambient to flight total pressure. On turbojet engines 
operating at low flight speeds, where np < ns, this term may be 
negative, indicating that a gain in thrust is made by passing part 
of the air through a pump and discharging it with an improved 
propulsive efficiency. This is equivalent to partially converting 
the turbojet engine to a turbofan and obtaining greater thrust at 
lower flight speeds because of the greater propulsive efficiency of 
the turbofan. 


Pressure Loss: 


A 
q 


This term represents the suction-drag coefficient associated with 
pumping the air from local static pressure on the wing up to 
ambient pressure. 


Duct Loss: 


Mo (1 — mp) 


This term represents the suction-drag coefficient associated with 
pressure loss in the duct between the wing slot and the pump. 

Curves of suction-drag coefficient versus flight Mach number 
are shown in Fig. 1 for both turboprop and_ turbo-jet 
engines. In each case the breakdown of the suction-drag co- 
efficient into the three components described above is shown. 
In these calculations, Ap/q is taken equal to 1.0 (corresponding 
to a fairly high lift coefficient), and 7p is taken to be 0.85. Both 
values are assumed to be independent of flight Mach number, 
Eq. (10) indicates that the suction-drag coefficient associated 
with pressure loss can be reduced in proportion to Ap/q, and that 
associated with duct loss can be reduced in proportion to 1 — np. 
The pump efficiency is taken as 0.8. 

The propulsive efficiency of the turboprop engines is assumed 
to be equal to the suction-pump efficiency, while that of the 
turbojet engines is given as the following function of flight Mach 
number: 

Flight Jet Propulsive 
Mach No. Efficiency 
0.3 0.33 
0.4 0.41 


0.5 0.48 
0.6 0.54 


The following conclusions can be drawn from Fig. 1. (1) Tur- 
boprop engines suffer a greater loss due to pumping than turbo- 
jet engines do. This is particularly true at low flight speeds. 
This does not necessarily imply that jet engines are superior to 
turboprop engines for boundary-layer-control applications; it 
does, however, indicate that the performance of turboprop- 
powered airplanes using boundary-layer control is more sensitive 
to the efficiency of the duct and pumping system than is the 
case for jet-powered airplanes. (2) At low-flight speeds, duct 
loss tends to be the largest contributor to the suction-drag 
coefficient. As flight speed is increased, pressure and momentum 
losses become more important. On low-speed turboprop air- 
planes particularly, care should be used in the design of the 
ducting. 

REFERENCE 
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